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MCDU-8  -  A  COMPUTER  CODE  FOR  THE 


NUMERICAL;  CALCULATION  OF  ONE -DIMENSIONAL  BLAST  WAVE  PROBLEMS 


INTRODUCTION 

The  propagation  of  blast  waves  in  an  inviscid  fluid  such  as  air  or 
water  has  always  been  of  interest ,  and  numerous  attempts  have  been  made 
to  obtain  its  solutions.  The  vast  majority  of  the  methods  used  required 
that  either  the  disturbances  be  weak,  or  that  the  explosions  obey  a  given 
similarity  constraint  which  is  appropriate  for  "point  source"  explosions 
only.  The  only  methods  which  can  give  complete  solutions  with  prescribed 
initial  and  boundary  conditions  are  numerical  methods  with  the  aid  of 
digital  computers. 

Among  the  numerical  calculations,  finite  difference  techniques  and 
the  method  of  characteristics  have  been  widely  adopted.  For  the  blast 
wave  problem,  shock  wave  propagation  is  one  of  the  most  important  features. 
Therefore,  the  method  of  characteristics  which  allows  shocks  to  be  traced 
exactly  is  inherently  more  accurate  than  finite  difference  methods  [1] . 

Although  there  are  many  research  papers  which  use  the  method  of 
characteristics  to  solve  blast  wave  problems,  most  of  them  have  made  re¬ 
strictive  approximations  for  simplicity.  Chou  and  Huang  [21  use  a  con¬ 
stant  time  scheme  in  conjunction  with  the  method  of  characteristics  to 
solve  a  blast  wave  problem  resulting  from  the  sudden  release  of  a  highly 
compressed  air  sphere.  Their  computer  code,  MCDU-7,  incorporates  a  strong 
shock  approximation. 

In  this  report,  the  numerical  method  and  computer  code  of  [2]  are 
modified  to  accept  any  equation  of  state  in  functional  form  involving 
pressure,  density,  and  specific  internal  energy.  A  technique  for  handling 
the  reflection  of  the  inward  traveling  shock  from  the  center  of  the  sphere 
is  also  included. 

In  the  first  section,  the  governing  equations  and  their  corresponding 
characteristic  equations  are  presented  followed  by  the  shock  equations. 

In  the  second  section  the  general  numerical  procedures  as  well  as 
details  concerning  the  calculation  of  certain  particular  points  are  describ¬ 
ed. 


In  the  third  section,  singularities  which  are  inherent  to  the  blast 
wave  problem  are  described  first.  The  solution  for  these  singularities 
then  follows. 
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In  the  fourth  and  fifth  sections  two  example  problems  are  solved 
and  compared  to  the  solutions  obtained  from  existing  computer  codes  and 
to  experimental  results.  The  first  problem  is  an  expansion  of  a  high 
pressure  sphere  with  an  ideal  gas  as  the  medium.  The  results  are  compared 
to  those  obtained  from  the  characteristic  code  MCPIJ-7,  which  is  restricted 
to  an  ideal  gas  medium.  The  second  calculation  is  for  the  explosion  of  a 
spherical  charge  50/50  Pentolite  (SO®  PFTN- 50"  TNT).  The  solution  is 
compared  to  those  obtained  by;  a.  the  Brinkley-Kirkwood  theory  [5]  , 
b.  another  computer  code  [4]  ,  and  c.  experimental  data  [5]  . 

In  the  sixth  section  some  conclusions  which  are  drawn  with  regard  to 
these  comparisons  are  presented. 

Appendices  I  and  II  contain  an  input-output  description  and  code 
listing  respectively. 


2 


I.  GOVERNING  EQUATIONS 


I- 

r 


The  governing  equations  for  one -dimensional  unsteady  motion  of  an 
inviscid  fluid  are: 


conservation  of  mass, 


3p  .  . 

3t  +  P 

—  +  u  —  +  (N-l)  p  -  =  0 

3r  ar  1  1  p  r 

(1) 

conservation  of  momentum, 

p 

au  .  au  A  ap  n 

3t  pu  F  Sr  =  n 

(2) 

and  conservation  of  energy 

3E  A 

at  u 

ar  p2l  at  u  ar 

(3) 

where  r  is  the  F.ulerian  space  coordinate;  t  is  time;  u  is  the  particle 
velocity;  p  is  the  density;  p  is  the  pressure;  and  F.  is  the  specific  in¬ 
ternal  energy.  In  these  equations,  N  is  a  constant,  with  values  1,  2,  and 
3,  corresponding  to  plane,  cylindrical,  and  spherical  waves,  respectively. 
Since  most  equations  of  state  are  given  as  a  relation  between  p,  p  and  E, 
we  shall  use  this  form  for  the  equation  of  state  in  our  calculations. 

Equations  (1),  (2),  and  (3)  are  a  set  of  first  order,  nonlinear, 
hyperbolic  partial  differential  equations.  Their  characteristic  directions 
and  equations  are 


along 


dr 

3t  = 

u  + 

c  & 

pC 

+  du 

+  (N-l) 

UC 

r 

dt  =  0 

(4) 

along 

dr 

at = 

u  - 

C  & 

PC 

-  du 

-  (N-l) 

uc 

r 

dt  =  0 

(5) 

along 

dr 

at  = 

c 

dE 

D 

'  p2 

dP  =  n 

(6) 

where 

the 

quantity 

c  is 

defined 

as 

c2  = 

(32.) 

F  +  E2  « 
t. 

^ap. 

(7) 

When  shocks  appear  in  the  flow  field,  the  material  properties  p,  E,  p 
and  u  are  discontinuous  and  the  above  mentioned  partial  differential  equations 
are  no  longer  adequate  to  describe  the  motion.  A  new  set  of  equations  is 
necessary  to  govern  the  propagation  of  these  shock  waves.  These  equations 
are 
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p2  ^'u2^  "  P1 
P2-Pi  =  Pi  GJ-Uj)  (u2_u1) 

VEi  *  I  ^v2)  (  ^  )  -  0 

where  U  is  the  shock  velocity,  and  the  subscripts  1  and  2  refer  to  the  states 
ahead  of  and  behind  the  shock  front,  respectively. 

II.  METHOD  OF  CALCULATION 

The  numerical  scheme  used  in  this  code  is  the  constant  time  scheme  uti¬ 
lizing  the  method  of  characteristics  introduced  by  Hartree  [6] .  This  scheme 
has  been  applied  by  Huang  and  Chou  [2]  for  the  calculation  of  expanding  high 
pressure  spheres  and  by  Chen  and  Chou  [1]  for  the  calculation  of  wave  propa¬ 
gation  due  to  intensive  in-depth  energy  deposition  in  a  two-layered  plate. 

They  showed  that  this  scheme  is  accurate,  and  can  be  easily  adapted  to 
computer  calculations. 

Although  the  governing  equations  used  in  the  present  code  are  different 
from  those  used  in  [1] ,  the  calculation  procedures  for  the  initiation  of  a 
second  shock  and  determining  the  properties  at  a  regular  point  in  the 
physical  plane  are  the  same.  The  starting  singularity  is  treated  "exactly" 
by  first  solving  for  the  properties  across  the  rarefaction  wave  and  then 
solving  for  the  properties  across  the  shock  wave  while  matching  the  pressure 
and  particle  velocity  at  the  interface  (contact  line) .  The  procedures  used 
are  quite  similar  to  those  found  in  [2],  thus  the  details  will  not  be  repeat¬ 
ed  here. 

The  arrangement  of  this  code  is  quite  different  from  previous  ones. 

This  code  consists  of  a  main  control  program  and  calculation  subroutines. 

Each  subroutine  is  designed  to  perform  a  specified  function.  These  subroutines 
may  be  classified  into  two  groups:  invariant  and  user  specified.  The  invar¬ 
iant  subroutines  need  never  be  changed  regardless  of  the  physical  problem 
or  materials  used.  For  example,  the  subroutines  for  calculating  the  proper¬ 
ties  at  different  types  of  points  in  the  physical  plane  are  common  for  all 
physical  problems  and  materials.  The  subroutines  that  must  be  changed  under 
different  situations  are  called  user  specified.  Each  of  the  major  subroutines 
and  its  function  will  be  spelled  out  as  follows: 

A.  Invariant  Subroutines 

General  point  subroutine  (Fig.  1):  Given  all  properties  at  the 
three  points,  II,  12  and  IS  along  a  constant  time  line,  this 
subroutine  calculates  all  properties  at  the  point  4  on  the 


(8) 

(9) 

(10) 


4 


next  time  line  by  using  the  three  characteristics  I,  II, 
and  III. 

4  4A  4R 


Figure  1 .  General  Point 


Figure  2.  Interface  Point 


Interface  (contact  line)  subroutine  (Fig.  2):  Given  all  proper¬ 
ties  at  the  interface  points  12,  13  and  the  neighboring  points 
II  and  14  all  on  a  constant  time  line,  this  subroutine 
calculates  the  new  interface  points  4A,  4B  on  the  next  time 
line,  matching  particle  velocity  and  pressure  at  both  points. 

Rarefaction  wave  subroutine  (Fig.  3):  At  the  initial  starting 
singularity  or  when  a  shock  reaches  a  free  surface  and  re¬ 
flects,  a  rarefaction  wave  occurs.  This  subroutine  calculates 
the  properties  across  the  rarefaction  wave  along  a  constant 
time  line  provided  that  the  pressure  behind  the  wave  and  all 
properties  in  front  of  the  wave  are  given.  A  parameter  must 
be  specified  to  control  the  number  of  subdivisions  that  each 
rarefaction  wave  will  be  broken  down  into.  By  increasing  the 
number  of  subdivisions  we  can  obtain  solutions  as  accurate  as 
we  wish. 


RARF  FACT  T  0\' 


Figure  3.  Shock  Reflection  From  a  Free  Surface 


Shock  equation  subroutine:  The  properties  behind  and  in  front  of 
a  shock  must  satisfy  the  shock  equations.  This  subroutine 
calculates  all  properties  behind  a  shock  given  all  properties 
in  front  of  the  wave  and  one  property  behind  the  wave. 

Shock  front  subroutine  (Fig.  4):  This  subroutine  calculates  the 
properties  behind  the  shock  front  at  point  4B,  provided  all 
properties  on  the  previous  constant  tune  line  and  one  condi¬ 
tion  behind  the  shock  front  is  known.  The  known  condition 
behind  the  shock  may  be  one  of  the  physical  variables  them¬ 
selves  or  may  be  given  in  the  form  of  a  characteristic  equa¬ 
tion.  For  the  latter  case,  the  characteristic  grid  is  shown 
in  Fig.  4A  and  Fig.  4B  for  right  and  left  traveling  shocks, 
respectively. 


Figure  4A.  Shock  Front  Point:  Figure  4R.  Shock  Front  Point: 

Right  Traveling  Shock.  Left  Traveling  Shock. 


Rarefaction  -  Shock  subroutine  (Fig.  5):  This  subroutine  solves 
for  the  initial  singularity  as  well  as  various  wave  inter¬ 
actions.  It  handles  two  possible  rases.  The  first  case 
happens  during  the  initial  explosion  (Fig.  5A)  or  during  the 
interaction  of  a  shock  wave  and  a  rarefaction  wave  (Fig.  5B) . 

It  consists  of  a  shock  traveling  to  the  right  and  a  rarefaction 
wave  traveling  to  the  left. 
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Figure  5A.  Starting  Singularity  Figure  SR.  Wave  Interaction  Point 

The  second  case  is  quite  similar  and  has  a  left  traveling 
shock  and  a  right  traveling  rarefaction  wave  as  shown  in 
Figs.  SC  and  5D.  An  iteration  procedure  is  used  to  match 


Figure  5C.  Starting  Singularity 


Figure  SO.  Wave  Interaction  Point 


& 
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& 
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the  particle  velocity  and  nressure  at  the  interface 
(contact  line). 

The  time  increment  calculation  subroutine:  This  subroutine 

calculates  the  time  increment  for  a  new  constant  time  line 
considering  numerical  stability  and  convergence.  We  have 
adopted  the  Courant -Friedrichs -Lewy  condition  as  the  stabi¬ 
lity  criterion  [7].  It  was  originally  derived  for  a  simple 
wave  equation;  however,  it  has  been  used  very  successfully 
for  more  complicated  sets  of  equations  [1],  [2],  and  [8]. 


n 
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The  point  arrangement  subroutine:  After  all  points  on  a  constant 
time  line  have  been  completely  calculated,  this  subroutine 
rearranges  the  points  before  proceeding  with  the  next  time 
line  calculations.  It  performs  the  following  steps: 

1.  Rearranges  the  order  of  points  with  respect  to 
their  position. 

2.  Automatically  adds  or  deletes  points  to  maintain 
a  relatively  constant  time  increment. 

3.  Automatically  deletes  those  points  whose  particle 
paths  cross  the  shock  waves. 

4.  Maintains  a  specified  number  of  points  in  front  of 
the  main  shock,  avoiding  the  calculation  of  points 
which  are  net  necessary. 

The  shock  reflection  subroutine:  This  subroutine  calculates  the 

properties  at  the  singularity  formed  as  a  inward  shock  re¬ 
flects  from  the  center  of  symmetry.  A  more  detailed  treat¬ 
ment  of  this  subroutine  will  be  given  in  the  next  section. 

The  center  point  subroutine:  This  subroutine  calculates  the  pro¬ 
perties  at  the  singularity  occurring  at  the  center  of 
symmetry  by  the  method  of  extrapolation. 

Initial  data  subroutine:  This  subroutine  assigns  all  properties 
to  all  points  along  the  first  constant  time  line  and  controls 
the  subroutines  used  to  solve  the  starting  singularity. 

R.  The  User  Specified  Subroutines 

Free  surface  subroutine:  This  subroutine  calculates  the  boundary 
point  of  a  physical  problem.  It  must  be  specified  for 
different  boundary  conditions. 

Non-dimensional  subroutine:  For  various  reasons,  it  is  beneficial 
to  non-dimensional ize  quantities  before  calculation.  This 
subroutine  must  be  adjusted  for  different  forms  of  non- 
dimensional  izat ion . 

Equation  of  state  subroutine:  This  group  of  subroutines  specifies 

the  equation  of  state  and  calculates  several  related  quantities. 
It  consists  of  six  subroutines;  EQSTCO,  EOSTEQ,  EQSTRO,  F.OSTPQ, 
EQSTPR,  and  EQSTPE  which  calculate  the  sound  speed,  c;  internal 
energy,  E;  density,  p;  pressure,  p;  and  the  derivatives  3p/3p 
and  3p/3E  respectively.  These  subroutines  must  all  be  changed 
for  different  equations  of  state.  Any  functional  or  tabular 
relation  among  the  density,  specific  internal  energy  and  pressure 
could  be  used  as  the  equation  of  state  in  this  code. 
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III.  SINGULARITIES 


In  either  the  expansion  of  a  spherical  compressed  pas  or  an  explosion, 
certain  mathematical  singularities  must  he  solved  before  the  general  numer¬ 
ical  calculations  can  begin.  If  we  start  calculation  from  the  instant  at 
which  the  detonation  wave  reaches  the  explosive  charge  surface  or  at  the 
moment  the  highly  compressed  gas  is  released,  a  discontinuity  of  properties 
exists  (the  so-called  starting  singularity).  We  solve  for  this  singularity 
by  using  regular  characteristic  methods  as  in  [2], 

At  the  center  of  symmetry,  the  particle  velocity  and  radius  are  both 
zero.  Tt  can  be  seen  from  the  equation  of  continuity  (1)  or  from  the 
characteristic  equations  (4)  and  (5)  that  the  term  u/r  becomes  uncertain  at 
this  point.  This  term  is  approximated  by  the  derivative  3u/3r  which  is  then 
extrapolated  from  the  neighboring  points. 

In  an  explosion  of  a  spherical  charge  there  is  a  shock  traveling  to¬ 
wards  the  center  of  the  sphere  in  addition  to  a  primary  strong  shock  propa¬ 
gating  outward.  This  inward  shock  is  usually  referred  to  as  the  second 
shock.  The  existence  of  this  second  shock  has  been  predicted  by  thcorv  [9]  , 
and  by  numerical  calculations  [10],  [11],  and  [12].  It  begins  as  a  very 
weak  compressive  discontinuity  which  builds  up  as  it  travels  toward  the 
center,  the  point  of  symmetry.  This  inward  traveling  shock  wave  will  have 
infinite  strength  (becomes  singular)  as  it  collapses  on  the  center. 

The  behavior  of  the  shock  near  the  singular  point  has  been  analytically 
studied  by  many  authors.  In  [13]  and  [14]  it  has  been  shorn  that  the  relations 
between  the  change  in  Mach  number,  M,  of  the  shock  wave  and  a  small  change 
in  the  cross-sectional  area,  A,  of  the  adjacent  particles  is  given  by  the 
formula 


.«■' —  a  n 

A  (PI-  -  11  K  (M) 

where  K(M)  is  a  slowly  varying  function  which  starts  at  0.5  for  a  weak 
shock,  M=l,  and  tends  to  0.394  as  M  ->  <»  (for  y=1.4).  Considering  a  point 
at  a  specified  distance  from  the  center  of  symmetry,  equation  (11)  shows 
that  the  Mach  number  will  be  the  same  regardless  of  whether  the  shock  is 
approaching  or  reflecting  from  the  center.  The  Mach  number  is  defined  to 
be  the  ratio  of  current  shock  speed  to  the  sound  speed  of  the  undisturbed 
medium.  The  sound  speed  of  the  undisturbed  medium  is  the  same  for  both 
reflecting  and  converging  shocks.  Therefore,  the  shock  velocity  of  both 
the  reflected  and  converging  shock  waves  at  any  arbitrarily  short  distance 
from  the  center  should  have  the  same  magnitude  but  be  opposite  in  sign. 
Using  this  conclusion,  we  next  present  a  brief  description  of  the  treat¬ 
ment  of  the  shock  reflection  as  used  in  this  code. 
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Let  us  assume  that  the  numerical  solutions  have  been  calculated  up 
to  a  time  ti  (see  Fig.  6),  the  time  just  before  the  shock  collapses  on 
the  center.  Examining  Fig.  6  we  can  see  that  the  converging  shock  inter¬ 
sects  time  line  t=ti,  at  r=ri .  Because  is  small,  we  assume  that  the 
shock  will  travel  the  short  distance  between  r}  and  the  center  with 
constant  speed.  Therefore,  the  location 


Figure  6.  Shock  Reflecting  From  Center  of 
Cylinder  or  Sphere 


and  the  velocity  of  the  reflected  shock  is  known  assuming  that  the  Mach 
number  of  the  reflected  shock  is  the  same  as  the  incident  shock.  The 
properties  across  the  reflected  shock  can  be  calculated  using  the  shock 
subroutine  and  assuming  that  the  shock  velocity  is  known.  After  solving 
for  this  singular  point,  the  solutions  of  all  other  points  can  be  calcu¬ 
lated  as  before. 
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IV.  A  SAMPLE  PROBLEM  INVOLVING  THE  SUDDEN  RE, LEASE  OF  A  IIICHLY  COMPRESSED 


AIR  SPHERE 


To  check  the  accuracy  of  this  code,  the  problem  solved  in  [2]  is 
solved  by  the  present  code  and  the  results  compared.  This  problem  in¬ 
volves  the  sudden  release  of  a  high  pressure  ideal  gas  sphere.  Before 
releasing,  the  properties  in  the  sphere  are  assumed  to  be  constant  with 
pressure  (P/Pa)  and  density  (p/pa)  ratios  with  respect  to  the  surround¬ 
ing  condition  of  100  and  1.16,  respectively.  The  specific  heat  ratio  is 
assumed  to  be  a  constant,  y=1.4,  for  both  media. 


This  same  problem  has  been  solved  by  the  present  code  and  the  results 
have  been  compared  to  the  ones  obtained  from  MCDU-7,  [2].  In  Fig.  7,  a 
comparison  of  the  physical  plane  produced  by  both  codes  is  presented. 

The  coordinates  x  and  X  are  the  dimensionless  coordinates  for  the  time 
and  radial  distance  from  the  center,  respectively. 

C  t 

t  =  -2—  and  x  =  — 

e  e 

where  Ca  is  the  constant  speed  of  sound  outside  the  wave  cone; 
c  =  [Fp/(10xPa)] J/3  is  length  expressing  energy  and  pressure  scaling;  Ep 
is  initial  total  energy  released;  and  Pa  is  the  constant  pressure  outside 
the  wave  zone.  For  this  sample  problem  with  the  initial  radius  of  the 
compressed  gas  of  1ft  (0.3048  m),  the  above  mentioned  quantities  are 


Ca  =  1116.7  ft/sec  (340.37  m/sec) 

Pa  =  14.7  psi  (1.013  x  105  M/m2) 
pa  =  0.07652  lb/ ft3  (1.226  kg/m3) 

Ep  =  2.19  x  106  ft-lb  (2.97  x  106  J) 


The  higher  shock  velocities  produced  by  MCDU-7  at  earlier  times  de¬ 
monstrates  the  effect  of  the  strong  shock  assumption  used  in  this  code. 
The  large  discrepancy  between  the  results  of  both  codes  for  the  inward 
traveling  shock  is  attributed  to  this  strong  shock  assumption.  Initi¬ 
ally,  the  inward  traveling  shock  is  much  weaker  than  the  outward  travel¬ 
ing  shock;  therefore,  the  strong  shock  assumption  will  result  in  larger 
errors  for  the  inward  shock  as  can  be  easily  seen  in  Fig.  8.  The 
solid  line  represents  the  results  calculated  by  the  current  code  and  the 
dotted  line  those  from  MCDU-7.  The  coordinate  tt  is  the  dimensionless 
pressure  P/Pa.  It  can  be  seen  that  the  strength  of  the  inward  shock,  S2 , 


INTERFACE 


Compressed  Air  Sohere. 


Figure  8  The  Pressure  Profiles  for  the  expansion  of  a 
Highly  Compressed  Air  Sphere. 


is  building  while  the  strength  of  the  main  shock,  Si  ,  is  decaying  as  time 
increases.  Between  the  two  shocks,  there  is  an  interface  where  a  discon¬ 
tinuity  in  the  slope  of  the  curve  occurs,  dust  after  r^n.075,  the  inward 
shock  reaches  the  center,  reflects,  and  propagates  outward  as  shown  in 
the  curve  for  t=0.086.  This  figure  also  shows  that  at  times  prior  to  reach 
ing  the  center,  the  wave  predicted  by  MCDU-7  is  faster  than  that  predicted 
by  the  present  code.  Again,  explanation  can  be  traced  to  the  fact  that  the 
initial  pressure  ratio  of  100  is  not  high  enough  to  justify  using  the 
strong  shock  approximation. 

Because  MCDU-7  cannot  handle  the  singularity  at  the  point  where  the  in 
ward  shock  reaches  the  center,  there  is  no  comparison  of  results  after  this 
t  ime . 


V.  A  SAMPLE  PRPBLfM  INVOLVING  A  BLAST  WAIT.  RFSirmC  FROM  THT  DETONATION 

OF  A  PHNTOLITF,  CHARGE 


The  second  example  problem  calculated  by  this  code  is  a  blast  wave  in 
air  produced  by  the  detonation  of  a  spherical  charge  of  Pentolite.  The 
blast  calculations  are  started  at  the  instant  the  detonation  wave  reaches 
the  surface  of  the  spherical  charge.  We  assume  that  the  resulting  gaseous 
products  of  the  detonation  have  reached  a  fixed  composition  as  we  start  the 
calculation. 


The  Abel  equation  of  state  [14]  is  used  for  the  explosive  products 


p  =  IJLX_£ 


a  p 


where  T  is  temperature,  R  is  the  universal  gas  constant,  N  is  the  number  of 
moles  of  gas  per  unit  mass  and  a  is  the  "co-volume"  of  the  gases.  For  this 
problem,  the  gaseous  products  of  the  solid  explosive  will  be  characterized 
by  an  ideal  gas  equation  of  state  (eq.  12  with  a=0) .  The  equation  of  state 
can  be  written  in  the  familiar  form  p  =  Hp(y-l).  The  constant  specific 
heat  ratio  of  the  explosive  gas,  ylt  is  calculated  from  the  properties  at 
the  detonation  wave  front  at  the  instant  it  reaches  the  charge  surface.  The 
value  of  Yj  used  in  this  problem  is  2.485.  The  medium  surrounding  the 
explosion  is  assuned  to  be  air  obeying  an  ideal  gas  equation  of  state  with  a 
constant  specific  heat  ratio,  y2 ,  of  1.4. 

The  data  concerning  the  conditions  when  the  detonation  wave  front  reaches 
the  charge  surface  are  obtained  from  [5].  The  properties  of  the  surrounding  c 
air  are  taken  as  standard  condition  at  sea  level,  i.e.,  pa  =  1  atm  (1.013  x  ID 
N/m2),  pa  =  1.293  x  10  3  gm/an3  (1.293  kg/m3). 
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For  convenience,  we  non -dimension!. ize  all  variables  before  calcula¬ 
tion.  The  dimensionless  variables  used  for  this  example  are  listed  as 
follows: 


p 


P 


p  C  2* 
M0  0 


F. 


f 


All  variables  with  a  "bar"  on  the  top  represent  dimensionless  quantities. 
The  reference  quantities  used  are  p0  =  2.58841  x  105  atm  (2.6227  x  1010 
N/m2) ,  and  p0  =  1  gm/cm3  (1.0  x  103  kg  /m3),  F,0  =  4.217  x  103  cal/gm 
(1.767  x  107  J/kg)  and  cq  =  8.0726  x  103  m/sec.  The  quantity  r0  repre¬ 
sents  the  charge  radius. 


The  wave  front  is  shown  in  a  physical  plane  plot  in  Fig.  9.  The 
coordinates  are  dimensionless  time,  t,  and  radius,  r.  A  second  shock  with 
zero  initial  strength  is  initiated  at  the  tail  of  the  rarefaction  wave 
very  early  in  the  calculations.  Although  the  second  shock  grows  in  strength 
and  propagates  inward  with  respect  to  the  explosive  gas,  the  large  particle 
velocity  of  the  gas  causes  the  shock  to  propagate  away  from  the  center  when 
viewed  with  respect  to  a  fixed  coordinate  system.  Consequently,  both  shocks 
in  Fig.  9  are  propagating  outward. 

The  relation  between  dimensionless  velocity  and  the  radius  is  presented 
in  Fig.  10.  It  can  be  seen  that  during  the  early  stages,  the  peak  value  of 
the  particle  velocity  is  high  and  concentrated  within  a  very  narrow  region. 
This  explains  why  in  the  early  stages  the  kinetic  energy  is  a  small  fraction 
of  the  total  energy.  As  the  wave  propagates,  the  contribution  of  the 
kinetic  energy  increases  and  the  potential  energy  or  internal  energy  de¬ 
creases.  The  two  discontinuities  in  the  velocity  profile  in  Fig.  10 
show  the  location  of  the  two  shocks.  It  also  can  be  seen  that  the  second 
shock  propagates  away  from  the  main  shock  front  as  time  increases. 
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The  shock  waves  in  the  Pentolite  Blast  Wave  Problem. 


Figure  11  shows  the  relation  hetw  .n  dimensionless  density  and  radius. 

At  any  instant  there  are  three  discontinuities  in  the  density  curves.  Two 
of  them  t.re  discontinuities  across  the  shock  fronts,  the  third  is  a  dis¬ 
continuity  across  the  interface. 

A  plot  of  the  dimensionless  pressure  distribution  with  respect  to  the 
radius  is  shown  in  Fig.  12.  (For  this  figure  and  Fig.  13,  pressure  has 
been  non-dimensionalized  with  respect  to  atmospheric  conditions,  i.e., 

P=P/Pa)  •  Initially,  when  the  detonation  wave  reaches  the  charge  surface, 
the  pressure  is  very  high;  approximately  a  quarter  million  atmospheres. 
Immediately,  a  shock  wave  forms  followed  by  a  rarefaction  wave.  The  pressure 
at  the  shock  front  drops  to  720  atmospheres.  In  a  very  short  time,  a 
compression  zone  appears  at  the  tail  of  the  rarefaction  wave  and  another 
shock  forms  (the  so-called  second  shock  or  inward  shock) .  Although  the 
strength  of  the  second  shock  is  growing  fast  and  it  is  traveling  to  the 
left  with  respect  to  the  particle  velocity  in  front  of  it,  the  absolute 
velocity  of  the  shock  carried  hv  the  explosive  gas  still  propagates  outward. 
As  time  goes  on,  the  second  shock  becomes  stronger  while  the  back  pressure 
becomes  lower  and  it  starts  propagating  away  from  the  main  shock  front. 
Fventually,  the  inward  shock  will  turn  toward  the  center  in  the  physical 
plane. 

Finally,  we  compared  our  results  to  those  obtained  from  an  existing 
code  [4],  Kirkwood- Brinkley  theory  [3]  and  experimental  data  [5].  From 
Fig.  13,  it  is  seen  that  for  early  times  our  code  gives  more  favorable  re¬ 
sults  in  comparison  to  the  experimental  data  than  others.  For  the  longer 
time  solution,  our  results  do  not  compare  favorably  with  the  experimental 
data.  We  obtained  a  pressure  ratio  across  the  main  front  shock  which  is 
higher  than  the  experimental  data.  This  may  be  the  results  of  using  a 
constant  specific  heat  ratio  in  our  calculations. 

The  computer  code  MCDIJ-8  has  been  run  on  both  the  IBM  360/7S  and 
the  Burrough  R5500  computers.  The  first  sample  problem,  took  approximately 
20  minutes  on  the  IBM  computer  for  all  results  shown  in  Figs.  7  and  8  with 
98  points  on  the  first  tune  line.  The  second  sample  problem,  took  80 
minutes  on  IBM  and  400  minutes  on  Burrouch  for  all  data  shown  in  Figs.  9  - 
13.  We  assigned  368  points  on  the  first  time  line  for  this  problem. 
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Fiqure  11  Density  Profiles  for  the  Pentolite  Blast 
Wave  Problem. 


Figure  13  The  comparisons  of  pressure  ratio  at  the  main 
shock  front  with  the  results  obtained  by  other 
existing  code,  Erinkley-Kirkwood  theory,  and 
experimental  data. 
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VI.  SlftMARY  AND  CONailSIONS 


A  one -dimensional  computer  code,  MCDU-8  has  been  developed  to  study 
the  problem  of  a  plane,  spherical  or  cylindrical  blast  wave  traveling 
through  an  inviscid  fluid.  This  program  uses  a  constant  time  scheme  in 
conjunction  with  the  method  of  characteristics  to  solve  for  the  flow  field 
in  regions  where  the  properties  behave  continuously  and  uses  the  Rankine- 
Hugoniot  relations  to  treat  shock  waves.  To  demonstrate  the  capabilities 
of  the  code,  two  sample  problems  have  been  calculated. 

The  first  sample  problem  treats  the  rapid  expansion  of  a  highly 
compressed  (100  atmospheres')  air  sphere.  The  calculations  begin  when  the 
air  is  released  and  extend  past  the  point  where  the  secondary  shock  wave 
reflects  from  the  center  of  the  sphere.  The  results  of  this  calculation 
are  compared  to  those  of  a  similar  characteristic  code,  MCDU-7  [2]  which 
utilizes  a  strong  shock  approximation  instead  of  the  more  exact  Rankinc- 
Hugoniot  relations.  This  comparison  shows  that  the  error  in  peak  pressure 
introduced  by  using  the  strong  shock  approximation  is  between  5  and  10 
percent  at  a  distance  from  tire  center  of  about  2.5  times  the  original 
radius  of  the  compressed  gas.  Due  to  the  limitations  of  MCDU-7,  we  were 
not  able  to  compare  the  two  codes  after  the  secondary  shock  reached  the 
center  of  the  sphere. 

The  second  sample  problem  treats  the  flow  field  produced  hv  the 
detonation  of  a  spherical  charge  of  Pentolite.  Like  the  previous  problem, 
the  calculations  begin  when  the  detonation  wave  reaches  the  surface  of 
the  explosive.  The  results  of  this  calculation  are  compared  to  those 
obtained  from  Kirkwood- Brink  ley  theory  [3],  experimental  data  [51,  and  a 
computer  code  developed  at  BRL  [4],  this  comparison  shows  that  during  the 
early  stages  of  computation,  MCDU-8  produces  results  that  are  in  closer 
agreement  to  the  experimental  data  than  the  two  other  means  of  calculation. 
It  has  been  concluded  that  if  one  is  only  interested  in  the  main  shock 
front,  then  the  Kirkwood- Brinkley  theoiy  is  adequate.  However,  if  detailed 
infonnation  concerning  the  entire  flow  field  is  of  interest  then,  the 
present  code  will  give  a  more  complete  analysis  than  other  available 
methods . 

It  is  hoped  that  in  the  future  we  will  be  able  to  extend  MCDU-8  to 
include  the  actual  detonation  calculations  thus  handling  the  complete 
explosion  problem. 
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APPENDIX  I  -  COMPUTER  CODE  DESCRIPTION 


A.  GENERAL  DESCRIPTION  OF  MCEIJ-8 

MCDU-8  is  written  in  single  precision  FORTRAN  IV.  It  is  designed  to 
numerically  solve  the  equations  governing  a  spherical  blast  wave,  by  using 
a  constant  time  step  iteration  scheme  in  conjunction  with  the  method  of 
characteristics . 

The  blast  wave  problem  consists  of  a  sphere  of  highly  compressed  gas 
(designated  region  one)  surrounded  by  another  gaseous  area  (designated  re¬ 
gion  two)  which  is  relatively  lower  in  pressure. 

MCDU-8  input  data  begins  with  a  card,  termed  the  'option'  card,  on 
which  the  user  selects  the  job  options  applicable  to  the  type  of  problem 
the  user  wishes  to  run.  MCDU-8  is  best  utilized  by  running  a  problem  to 
completion  in  a  series  of  small  computer  runs  rather  than  all  at  once.  The 
option  card  provides  an  easy  method  for  doing  this.  After  reviewing  the 
output  from  any  particular  computer  run,  the  user  has  three  options: 

1)  The  problem  may  be  carried  out  to  a  larger  time  using 
the  same  time  step  as  was  used  in  the  previous  run; 

2)  the  problem  may  be  carried  out  to  a  larger  time  using 
a  different  time  step  or, 

3)  the  same  run  may  be  repeated  using  a  different  time  step. 

The  first  two  options  give  the  user  control  over  the  rapidity  with  which 
MCDU-8  calculates  a  solution  while  the  third  option  allows  the  user  to  improve 
the  accuracy  of  any  particular  run. 

B.  INSTALLATION  DEPENDENT  FEATURES 

MCDU-8  utilizes  two  data  files  which  must  be  made  available  to  the  pro¬ 
gram.  The  opening  of  data  files  and  the  devices  on  which  they  are  stored,  for 
example,  tape  and  disk,  is  a  function  of  the  job  control  language  and  the  faci¬ 
lities  available  at  any  particular  installation.  The  user  should  insure  that 
the  two  files  have  the  following  characteristics: 

1)  The  files  must  have  the  unit  numbers  one  (1)  and  two  (2); 

2)  they  must  have  a  physical  record  size  of  at  least  five  (5)  words; 

3)  a  minimum  record  length  of  2010  records  (10,050  words)  and, 

4)  all  I/O  is  unformated  and  performed  serially  (random  access  is 
NOT  used) . 

C.  OPTION  CARD 

MCDU-8  always  requires  at  least  one  card  of  input,  termed  the  'option' 
card,  on  which  the  user  specifies  what  actions  the  program  is  to  take  in 
solving  a  problem.  This  card  is  always  the  first  card  of  the  input  deck. 

The  six  variables  initialized  by  this  card  are  listed  in  Table  1  for  refer¬ 
ence.  Examples  of  the  option  card's  uses  are  given  in  later  sections. 
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TARLH  1  OPTION  CARD  VARIABLES 


VARIABLE  COLUMNS  EDRMAT  DESCRIPTION 

ISTART  1-2  12  =  -1  (Continue  the  problem  with  card 

input) . 

=  0  (Start  a  new  problem  with  card 

input) . 

=  1  (Continue  the  problem  with  file 

one  input) . 

=  2  (Rerun  previous  run  with  file 
two  input) . 


I  PI  INCH 

3 

11 

=  0  (No  punched  output) . 

=  1  (Punch  out  the  last  calculated 
time  line.  These  cards  arc  used 
with  ISTART  =-l  to  continue  a 
problem) . 

I  DUMP 

4-6 

13 

(Calculated  time  lines  are  periodic¬ 
ally  stored  on  file  one  for  safe  keep 
ing.  13  specifies  how  many  lines  are 
to  be  calculated  before  a  line  is 
dumped  onto  file  one). 

I1TT 

7 

11 

=  0  (No  new  time  step). 

=  1  (New  time  step  is  to  be  specific* 
in  columns  8  through  22) . 

DT2 

8-22 

E15.8 

(New  time  step  to  be  used  if  TPMl. 

TMAX 

23-37 

E15.8 

(Problem  time  to  which  MCDII-8  is  to 
calculate  a  solution) . 

D.  DEFINING  A  NEW  PROBLEM 


Figure  14  illustrates  the  cards  needed  to  define  a  new  problem. 

Besides  the  option  card,  six  additional  input  cards  are  required. 

The  first  card  in  Fig.  14  is  the  option  card.  A  new  problem  is 
signaled  by  setting  ISTART=0.  IPUNCH  has  been  set  equal  to  one.  Upon 
completion  of  a  computer  run,  the  final  time  line  will  be  punched  out. 

The  user  need  only  include  an  option  card  at  the  beginning  of  these 
punched  cards  with  ISTART=-1  in  order  to  continue  a  run.  IDUMP  has  a 
value  of  three,  thus  every  third  time  line  calculated  will  be  dimped 
onto  file  one.  When  an  option  card  is  encountered  with  ISTART=1,  the 
last  time  line  dumped  onto  file  one  will  be  read  and  used  to  continue  a 
problem. 

When  a  problem  is  first  defined,  a  singularity  exists  at  the  inter¬ 
face  between  regions  one  and  two  because  of  the  difference  in  pressure. 
MCDIJ  8  has  the  ability  to  pick  its  own  time  step.  The  user  may  specify 
a  time  step  by  setting  HTT=1  and  placing  the  time  step  in  columns  8 
through  22.  Accidently  specifying  too  large  a  time  step  may  cause 
serious  difficulties  with  the  program  logic.  As  in  example  2  ,  it  is 
usually  good  policy  to  let  MCDU-8  calculate  the  time  step  and  to  run 
the  problem  out  several  time  lines  until  the  rarefaction  wave  occurring 
at  the  singularity  becomes  "smeared”  through  out  region  one.  This  is 
done  by  setting  IDT=0.  Notice  that  TMAX  has  been  set  equal  to  fifty 
microseconds . 

Card  two  contains  the  initial  time  step  to  be  used  to  calculate  the 
properties  about  the  singularity.  This  consists  of  a  shock  travelling 
into  the  region  two  and  a  rarefaction  wave  with  velocity  towards  the 
center  of  region  one.  For  good  results,  a  time  step  should  be  chosen  so 
that  after  the  calculation  of  the  initial  singularity  about  the  interface, 
the  head  of  the  rarefaction  wave  has  traveled  approximately  three  percent 
or  less  of  the  distance  from  the  interface  to  the  center  of  region  one. 

In  this  region,  the  shock  strength  is  assumed  constant  and  properties 
remain  constant  along  characteristic  directions.  Following  these  assump¬ 
tions,  the  code  assigns  mesh  points  at  locations  where  the  characteristics 
emulating  from  the  singularity  cross  the  first  time  line.  The  code  then 
uses  stability  criterion  to  determine  the  time  increment  to  the  next  time 
line.  Format  is  E15.8. 

Card  three  supplies  the  specific  heat  ratios  (GAMMA)  oF  the  two 
regions.  In  the  example  given,  these  values  are  both  1.4.  Format  is 
2E15.8. 

Card  four  lists  the  dimensions  of  the  data  supplied  to  the  program. 
This  is  only  a  programming  convenience  for  the  user  and  will  label  the 
output  for  documentation  purposes.  If  this  card  is  left  blank,  no  labeling 
will  result. (see  Table  2). 
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TABLE  2  -  CARD  4  -  OUTPUT  LABEL 


QUANTITY 

EXAMPLE 

COLUMNS 

Time 

Seconds 

1-10 

Particle  Velocity 

Ft/Sec 

11-20 

Sound  Speed 

Ft/Sec 

21-30 

Pressure 

Lb/Ft**2 

31-40 

Mass  Density 

Slug/Ft* *3 

41-50 

Specific  Energy 

Ft -Lb/Slug 

51-60 

Length 

Ft 

61-70 

Card  five  defines  the  physical  properties  or  region  one:  pressure 
CPU) ,  particle  velocity  (UI1)  and  density  (RHI1).  In  the  example  of 
Fig.  14,  pressure  is  equal  to  2.1168  105  lb/ft2,  particle  velocity  is 
zero,  and  density  is  2.7739  10-3  slugs/ft3.  Format  is  3E15.8. 

Card  six  contains  the  physical  properties  of  region  two.  Pressure 
(PI2)  in  Fig.  14  is  2.1168  103  lb/ft2 ,  particle  velocity  (IJI2)  is  zero, 
and  density  (RHI2)  is  2.3913  10“ 3  slugs/ft3.  Format  is  3E15.8. 

Card  seven  contains  three  values.  The  first  value  supplied  on  this 
card,  (IN),  specifies  the  number  of  subdivisions  the  initial  rarefaction 
wave  is  to  be  divided  into  for  calculation  purposes.  In  the  example  of 
Fig.  14,  IN  is  set  equal  to  nine.  Ideally,  the  distribution  of  pressures 
from  the  first  point  in  the  rarefaction  wave  to  the  last  point  should  be 
smooth  and  free  from  large  jumps.  Only  experience  can  determine  what 
value  of  IN  will  give  a  reasonable  solution  to  a  problem.  Theoretically, 
we  can  make  the  solution  as  accurate  as  we  wish  by  choosing  large  values 
of  IN.  Practically,  IN  should  never  exceed  twenty-nine  as  MCDU-8  can 
only  handle  up  to  this  many  subdivisions  without  enlarging  the  storage 
allocation.  Format  is  12. 

The  second  value  supplied  on  this  card,  (XZ) ,  is  the  initial  radius 
of  the  sphere.  Format  F.15,8. 

The  third  value  supplied  on  this  card  is  the  pressure  jump  toler¬ 
ance,  (PTOL) .  FTOL  specifies  what  fractional  percent  pressure  rise  per 
unit  distance  must  be  present  between  two  points  before  a  shock  wave  will 
be  initiated.  MCDU-8  has  the  capability  to  initiate  one  left  traveling 
shock  wave.  In  Fig.  14,  this  has  a  value  of  4  which  is  equivalent  to  a 
4000  percent  pressure  jump  per  foot.  This  value  appears  to  yield  satis¬ 
factory  results  for  this  problem.  This  value  may  need  to  be  adjusted  to 
suit  a  particular  problem.  Format  is  H15.8. 
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E.  RESTARTING  FROM  FILE  ONE 


The  IDUMP  variable  on  the  option  card  specifies  how  many  tine  lines 
arc  to  be  calculated  between  dumps  on  file  one.  If  IDUMP  is  set  equal 
to  zero,  no  time  lines  will  be  dumped.  Each  new  time  line  dumped  on  file 
one  replaces  the  previous  one.  If  after  a  run,  the  user  decides  to  carry 
calculations  out  to  a  greater  time,  the  user  need  only  to  read  in  an  op¬ 
tion  card  with  ISTART:T.  The  last  line  dumped  will  be  read  off  of  file 
one  and  calculations  will  proceed  from  there. 

By  setting  I0T=1,  the  user  may  also  change  the  time  step  by  supplying 
the  new  time  step  in  columns  eight  through  twenty- two  on  the  option  card. 

In  the  example  of  Fig.  15,  a  restart  fro  file  one  has  been  called  for 
(ISTART=1),  no  punched  output  is  requested  (IPUNCH=0) ,  every  time  line  will 
be  dumped  (IDUMP=1) ,  a  new  time  step  is  called  for  (IDT=1,  DT2=1.0  10'6 
seconds),  and  the  problem  will  be  calculated  out  to  a  time  of  100  micro¬ 
seconds  (TMAX:  100.0  10‘ 6  seconds). 

An  added  feature  of  the  file  one  restart  is  if  a  job  terminates  ab¬ 
normally  (such  as  running  out  of  computer  time)  and  the  program  was  not 
executing  I/O  with  file  one  (causing  parts  of  two  different  time  lines  to 
be  saved)  the  user  may  restart  the  problem  as  explained  above  and  MCDU-8 
will  proceed  from  the  point  of  termination. 

F.  RESTARTING  FROM  FILE  TOO 

Whenever  a  new  job  is  begun,  the  very  first  time  line  calculated  is 
dumped  onto  file  two.  A  user  may  repeat  the  same  run  over  as  many  times 
as  desired  by  simply  reading  in  an  option  card  with  ISTART=2.  (See  Fig.  16) 
by  changing  the  time  step  and  comparing  the  results  with  earlier  runs,  the 
convergence  of  the  problem  to  an  accurate  solution  can  be  checked.  In  the 
example  of  Fig.  16,  the  time  step  has  been  set  to  one-half  microsecond. 

Remember  that  whenever  a  new  run  is  started,  the  first  time  line  of 
the  run  replaces  the  first  time  line  of  the  old  run. 

G.  RESTARTING  FROM  PUNCHED  OUTPUT 

If  several  problems  must  be  run  on  MCDU-i  concurrently,  it  becomes 
impractical  to  provide  disk  or  tape  restart  capabilities  for  every  problem. 
MCDU-8  can  remember  time  lines  for  only  one  problem  at  a  time.  By  setting 
IPUNCH=1  on  the  option  card,  the  last  time  line  calculated  at  the  end  of  a 
run  will  be  punched  out.  One  need  only  place  an  option  card  with  ISTART=J 
at  the  beginning  of  this  punched  output  deck,  being  careful  to  preserve  the 
order  of  the  cards,  in  order  to  continue  calculations  from  the  end  of  the 
previous  run  of  that  problem.  All  of  the  other  features  of  the  option  card 
are  still  available  to  the  user  of  MCDU-8.  (See  Fig.  17). 
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II.  SAMP  LI;  OUTPUT 

Figures  18  thru  23  present  sample  output  for  MCDU-8.  The  information 
is  completely  labeled  and  should  he  self  explanatory. 


u;  1  J  U 

_»  J  T  T  «  - 

T  U  ^  —  ■* 
I  — 

►  r  ii  *tu 

*  E  5  W  —  *1 

—  *  *3  X  « 

3  1  i/>  1  I  </1 


»  I 
t  Jl 

<  ft  * 

*-3  0*-  * 

irt  1  3  3  X 

-3  ►- 


BBT  available  copy 


Figure  18.  Option  Card  Printout. 


Figure  20.  Initial  Singularity  Printout  Generated  by  A  New  Problem. 
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Figure  22.  Typical  Tine  Line  Printout  (Continued) . 
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C  -  MAIN  PROGRAM  FUR  BLA&T  WAVE  * 

FILF  l«BLASTA/ME?367»UNlT«0lSK»BLQCKlNG«30»R£CURD«b 
FiLF  2«BLASTR/M£?367*UNlT*QlSKiBLQCKlNG*30#REC0RD*b 

common /gain/  q(2#iooo) »x(2. iooo).u(2»iooo) .C!2» ioooj.rp!?. i 

1 OOO  )  *  E!?»1000).P!2»1000) 

COMMON/TIMUU/  DT  »  UU1»  UU2 »  I » XMU 

COMMON/NOCON/  IS1.  IS?*  I S  3  *  ISA*  INTI*  INT2*  IMAX 

COMMON  / I N I T /  PI1»UI1»RhI1*EI1.CI1*QI1»PI2»UI2»RHI?»EI?»CI? 


vol 

902 

903 

904 

905 

906 

907 


1*012 

Common  /  refl  /  tref.t 
commqn/otts/  ott 
common/no IM/TT » i tmax»xaz#Tmax 
COMMON/SHKl/  £P 
C0MM0N/GAM/GAMM(2) 
C0MM0n/C0N/TC0N»DT2*PTUL 
C OMMUN/NNEWW/K SMOCK »KTtLG* IDT 
FORMAT! 15. 3E15.8) 

roRMATUH  ,60!/."  ".7!"*")*6!"MCDU8»QLD 

** 


" ) »  7! "* 
"  ) »  7 !  "  * 


"  )  )  ) 
"  )  )  ) 


"*  IPElS.rt/ 
".1PE15.B/ 
".  lPtlb.B) 


PROBLEM 

FUHMATUH  .60!/."  ",7l  %").6!"MCDU8.NEW  PROBLEM 
FORMAT!  1H1 ."LAST  TIME  ^  I NE  WILL  BE  PUNCHED") 

Format (sei5»8) 

FORMAT! 1615) 

FORMAT!  1HQ»"NUN-DIMF.NS*0NAL  INPUT  DATA"/ 

!1H  #"!T)TIME  OF  THE  INITIAL  TIME  LINE 
11H  ."!OT)TImE  iNCHEMENf  To  THE  FIRST  TIME  LINE 
(1H  »"! TMAX ) M A  X  I  MUM  RUN  TIME 

908  FORMAT! 1H0»"GAMMA  FOR  REGION  ONE  ".1PE15.8/ 
llH  » "GAMMA  FUR  REGION  f wO  "#1PE15.8) 

909  FORMAT! 1  HO . " ! UU 1 ) LEF T  TRAVELING  SHOCK  VELOCITY  ".1PE15.8/ 
!1H  *"!UU2)RIflMT  TRAVELING  SHOCK  VELOCITY  ".1PE15.H) 

910  FqRMAT!T57*I4.T1»I5»E1^»8) 

911  FURMATUH  .T57.I4.Tl.Ib.E15. 8) 

913  FORMAT! 1H1 ."PROPERTIES  OF  THE  INITIAL  TIME  LINE"/ 

! 1  HO . "POINT  NO.  "  »7X."*",14X."U".14X."C"»l4X»"P".13X."RH"»l 
14X»"E") 


914  FORMAT! 1H  . 3 X . 1 4 . 3 X . 6 !  1  PE  1 5 . 8 ) ) 

9 ?7  FORMAT! 1H1 .6!"  MCOUH.TlME  LINE".I4)/ 

*1H  ."TIME  OF  THIS  LINE  ".lPElb.8/ 

iiH  ."time  increment  to  previous  time  line  ".ipeis.b/ 

!  1  H  ."LINE  REGION  PO I  NT” . 7X . "X" » 1 4X . "U" . 1 4X . "C" . 1 4 X . "P" » 1 4 X . 
1  "RH" »  13X."E") 

9?8  FORMAT! 1H1 ."LINE  REGION  PO INT" . 7x » "X" . 1 4X » "U" . 1 4 X  .  "C" . 1 4x . " 
IP".  13X»"KH"»14X»"£") 

9?9  F0RMAT11H  ."ThE  NEXT  T*0  POINTS  DEFINE  THE  LEFT  TRAVELING  S 
l HOCK" ) 

930  FflRMATUH  ."THE  NEXT  TWO  POINTS  DEFINE  THE  INTERFACE") 

931  FflRMATUH  ."THE  NEXT  TWO  POINTS  DEFINE  THE  RIGHT  TRAVELING 
1  SHOCK" ) 

932  FORMATUH  .  1 4 . 1 X .  1 6 . 1  X  »  1 5 . 6!  1PE  15 . 8 ) ) 

933  FORMAT! 1H1.601/.12!"  E^D  *  MCDUB" ) ) ) 
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914  FORMAT! 1H1#"PUNCHED  OUTPUT  OF  THE  LAST  TIME  LINE" ) 

915  FqRMAT<T57*I4#T1»3E15.») 

916  F0RHAT(T57#I4»T1»7I5) 

918  Formats  this  card  to  contain  istart#ipunch»idump»idt»ot.tm 
1AX") 

919  FoRMATCIH  »T57»l4»T2*3tl5,8) 

94O  FORMATUH  #  T57 # 1 4  #  T2 »  7  I  5  > 

94  1  FuRMAT(///lH  # tOX#"TIMt< T>"#2X#"*"#30X»"*"#2X#"INTFRFACE"/ 
(1H  #18X#"***"#28X#"***"» 

( 5( /# 1H  .19X*"*"»30X»"*")/ 

(1H  #19X#"*  REGION  l"#'52#"*  REGION  2"# 

C  4  (  /  »  l  H  #19X#"*"»30X#"*")/ 
llH  #19X#"*"#30X#"*"#29*#"*"/ 

(1H  #  19X  »63( •-  DISTANCE  FROM  THE  CENTER  POINT  (X)"/ 

ClH  #80X#"*"//// 

!1hO#"TITLE  ABBREVIATIONS"/ 
t 1  HO  # "X  -DISTANCE  FROM  ThE  CENTER  POINT"/ 

ClH  #"U  -PARTICLE  VELOCITY"/ 
llH  *"C  -SOUNO  SPEED"/ 

(  1H  »"P  -PRESSURE"/ 

ClH  #"KH-MATERIAL  DENSITY"/ 
l 1H  # "E  -SPECIFIC  ENERGY") 

942  FuRMATC 1  HO # "  (  UU2 ) R 1 GHT  SHOCK  VELOCITY  "#1PE15.8) 

943  F  Q  R  M  A  T  ( I2*I1*I3»I1»2£1^«8) 

944  FORMATUH  #"(UU1)LEFT  ^hOCK  VELOCITY"*  1PE1S»8) 

946  F ORMAT ( //// 1 H  ."(ISTARl)  "#I3/1H  # 

("(IPUNCH)  " *  I  3/ 1 H  » 
l " ( I  DUMP  )  ",13/1H  » 

("(IOT  )  " #  I  3/  1 H  # 

("(OT  )  ".1PE15.8/1H  » 

("(TMAX  )  ".1PE15.8) 

TrEF  *  10, 

XMU  *  3, 

IREF  «  l 
I  ■  1 
L  *  1 
K0UMP«0 

HEAD  94 3, 1  ST ART* IPUNCH*  I OUMP # IOT » DT2 » TM AX 
I F  C ISTART.LT.IJGOTO  1 
PRINT  902 
PRINT  941 

PRINT  946* I  START# IPUNC*#  IDUMP# IDT* DT2* TMAX 
CALC  KEaOO! ISTART) 

CALL  OUMP( 2 ) 

GO  TO  8054 
1  CONTINUF 

IF(  ISTART. LT.OGO  TO  2 

print  903 

PRINT  941 

print  946 » i st art » i punch »i dump »idt»dt2»tm ax 
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call  inidat 

KTELG-l 
KSHOCK**  1 

call  dump(2) 

GO  TO  80S3 

2  CONTINUE 

print  902 

print  941 

PRINT  946# IsTART  » I PUNC  R  # I  DUMP  # IDT  #0T2»TMAX 

HEAD  90S »T#DT#TC0N 

tmax»tmax/tcun 

IF(  IDT. EG.  1  )0T*DT2/TC0'>* 

R£ AO  90S.GAHMC 1 ) »GAMM( *) 

HEAD  90S  #  UU 1 #  UU2 
READ  906»lSl#IS2#lS3#IaA»IMAX 
READ  90fr.INTl»lNT2#KTE‘*G»KSHnCK 
print  9o7.t,ot#tmax 
print  908#gamm( i j #gammi?) 

Print  909#uui»uu? 

PRINT  913 

DO  3  J-l.IMAX 

REAU  90S.X(1.J).U(1,J)#C(1»J) 

REAU  90s»P(l»J)»RH(l,j)»ECl»J) 

0(1#  J)«0*0 

PRINT  91A»J.X(1»J)#U(1'J)#C(1»J)»P(1#J)#RH(1#J)»E(1»J) 

3  continue 

REAU  90 1 » I » PTOL 
Call  dumP(2) 
bosA  continue 
TnEn»ut 

5  CUNTINUE 

IFUSHUCK.EQ. 1 )G0  TO  80 
1FC UU1  .GT.O.OJGU  TO  80 

lF( IS1  #LE.3#UH.UT*UUl  +  *( 1 # IS1  )  #LT.X( 1 » 3) )GD  TU  8 

ho  continue 
call  ptarng 
go  to  6 
8  CONTINUE 
EP  rn  0#0l 
CALL  REELSK(l) 

IRE.F  a  2 
TREF  ■  T 
T  ■  T  ♦  OT 

CALL  SHFPT( ?# ISJ-2#  IS 3*1 »IS3# IS4. 1 S4^ 1 > IS«  +  2#UU?) 

GO  TO  12 

6  continue 
t.t+dt 
TnEw  ■  nT 


TeP  ■  T*.2l5 

i r (  k shock i rtt « 2 )  go  to  «s5oi 

IS2P2*? 

DO  6500  J* 1 » I  NT  1 *2 
C 
C 

c 

C  PhEGSURE  slope  test 

TEST*IP(1»J*1)-P(1»J))'CPC1»J)*IX(1»J*1)-XC1»J))) 

IF(  TEST.LT. PtUDGU  TO  6500 
c 
c 
c 

lF((U(i,J)-CU*J))»LE.<U(l*J  +  I)-C(l»JM)))GQ  TO  6500 
I  S  1 »  J 

PRINT  650?* I S 1 

6502  FORHATUH  ."A  SHOCK  w I *- L  BE  INSERTED  AT  POINT  NUMBER  "*I5) 

I  S2*  I S 1 ♦ l 

Call  smitcht ist  +  i »i»i *d 

InTi*intui 

InT2*INT2*1 

1 S 3*  I  S3* 1 

IS4*IS4*1 

Ep*.  1*(P(  1»ISU2)-P(  1MS1))/PC1»IS1) 

call  shokini i *2* i » isi  ) 

kshuck*? 

00  TO  6501 

6500  CONTINUE 

6501  CONTINUE 

I F  t KSHOCK  »EU . 1 )  00  TO  6504 
CALL  SHFPT( 1 » IS1-2# IS1“1 » ISi * IS2. IS2+1 . IS2+2.UU1 ) 
IF(P(2»IS2).LE.P(2»IS1 > )KSH0CK«1 
1 F ( KSHOCK • Eq • 2 )  GO  TO  H001 
PRINT  6503 

6503  FoRMATdH  ."DUE  TO  OE'CAY  IN  PRESSUKE*THE  INSERTED  SHOCK  HAS 
1  BEEN  RFMOVEU") 

call  switch(IS2M#-i*i»2> 

1nT1«INT1-1 
InT2«INT2*1 
I S  3* I  S3* 1 
I  54* I S4*l 
I S  1*4 
I  S2*4 

4001  CONTINUE 

IS1M2  •  ISI  -  2 
DO  10  K  *  2.  I S 1  M2 
10  CALL  UNPTll#  K* 1 »  K  *  K  +  1#  3) 

IFCUU1  ( GE •  0. )  Kfl  ■  ? 

1F(UU1  ,LT,  0. )  KQ  ■  3 
KQ  *  J 
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CALL  GNPT ( t (  IS1-2.  I S i -  1  *  IS1»  KQ) 

12  IF(IS2  ♦  1  .EQ«  INTI)  on  TO  21 
1 F  C  UU  1  . Gfc  •  U. )  KU  ■  a 
1 F ( UU 1  ,LT,  o.)  KU  ■  1 
KU  *  J 

1F(KSH0CK.EQ. 1 )  CALL  GNPT(t»lSl*l*ISl*ISl+l*3> 

Call  GNPT ( 1 #  1 S2  »  I S2+ 1  *  IS2+2*  KQ) 

15  CONTINUE 

IS2P2  ■  IS2  ♦  2 
6504  CONTINUE 

INT1M1  *  INTI  -  1 

IF(IS2P2  .GT.  INT1M1  )  GO  TO  21 

00  20  K  *  I S 2 P 2 »  I  NT  1 H  1 

?o  call  gnptci.  k-i«  k»  k  +  i*  3) 

?1  Call  INTFPT( I  NT  1  - 1 .  IN'1»  INT2#  INT2*1) 

1S3M2  »  Ii>J  •  2 
1 NT2P  1  >  I  NT 2  ♦  1 
IF  IInT?P1  ,GT.  IS3H?  )  GO  TO  32 
UO  30  K  ■  INT2P1*  IS3M2 

30  Call  gnpt  < ?»k-i »k»k*i* 3) 

IFIUU2  .  GE «  0.  )  KU  «  2 
IFIUU2  .LT .  0. )  KU  »  3 
KO  ■  i 

12  CALL  uNPT(2.  IS3-2#  IS-i-1*  ISi»  KU) 

IF  (  UU?  .GF.  0.  )  KO  *  3 
IF  (  UU?  .LT.  0.  )  KO  3  1 
CALL  ShFPT(?.IS3-2.ISJ-1#IS3»IS4»IS4*1#IS4+2»Uu2) 
Kg  *  3 

CALL  GNPT  (?.IS4.IS«M»IS4*2»KQ) 

IS4P2  *  IS4  ♦  2 
1MAXM1  m  I M A X  -  1 
00  40  K  *  IS4P2*  IHAXMl 
aO  CALL  GNPT(2.  K-l.  K#  K*l»  3) 

1 F ( T  .Eo.  THEf)  GU  TO  “2 
call  cEntptcw  i#  2.  3> 

42  CONTINUE 

Call  AD0PTS(2) 

I  *  I  ♦  1 

1)0  45  J  =  1,  I H A  X 

X  t  i  »  J  )  «  X  (  2  *  J  ) 

P{ 1#  J)  *  P ( 2 *  J) 

U ( 1 .  J)  *  U( 2#  J) 

HH( 1 .  J)  *  RH ( 2  »  J) 

£( 1#  J)  »  E ( 2 *  J ) 

C ( 1  *  j)  ■  C( 2#  J) 

U ( 2  *  J)  *  0. 

45  U(  1 #  J  )  a  Q ( 2 •  J ) 

bo«si  continue 
temp.ut 

I F ( I.EO.l )TEMP-T 
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KDUMP*K0UMP*1 

IF(KOUMP,NE.IOUMP)GG  TU  120 
CALL  DUMP(l) 

KOUHPaO 

i ?o  continue 

PRINT  927»I»I»I»I»ItI*T(TEMP 
KCUUNT*0 
Of)  950  J*1»IHAX 
kcount*kcqunT+i 
11*1 

IF( J«GE«  INT?) 1 1*2 

IF( J.EQ.IS1.ANU.KSHQCK*EQ<2  )  PRINT  929 

IF(  J.EO.INTDPRINT  930 
lF(J#Efl«IS3) PRINT  931 

PRINT  932»I»II.J#X(1#J).U(1*J)#C(1#J)»P(1*J)#RH(1#J)»£(1#J) 
lF(KCUUNTiNE.55)GU  TO  *50 
KC0UNT*0 
PRINT  928 

950  continue 

PRINT  9«2#UU2 
I F  C KSHOCK  #E9 » 1 )  GO  TO  945 
PRINT  944  »  UU 1 
94S  CONTINUE 

IF( T  *GT «TMAX )G0  TO  50 
IF( I »E0» 1 )G0  TO  8054 
GO  TO  5 

so  continue 

IF(  IPUNCR.Eu*0)G0  TO  9<>1 

9s5  continue 
i  i- i 

print  934 
1*1 

punch  93a 
PRINT  938 
I  *2 

PUNCH  935# I »T#OT»  TCQN 
print  939*i#t.ot.tcon 
1*3 

PUNCH  93S»I»GAMM(  1)»GAmM(2) 
print  939#i»gamhc i),gaMm<2) 

1*4 

PUNCH  935  » I  *  UU 1 1 UU2 
PRINT  939»I»UU1»UU2 
1*5 

PUNCH  936»I»IS1» IS2»IS3,IS4#IHAX 

print  940*1. isi»is2#is4#is4»ihax 
1*6 

PUNCH  936» I »INT1'INT2»KTELG» KSHOCK 

PRINT  940#I.INT1.INT2»RTELG# KSHOCK 
00  960  J* 1 # I H AX 
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I* I  ♦  1 

HUNCH  93S#  1  ,X(2, J) ,U( ?• J) »C< 2, J> 
PRINT  939»I»X(2»J)»U(2,J)»C(?»J) 

1  >  I  ♦  1 

PUNCH  9  35.I,P(2.J),RK(<?»J)«E(2,J) 
PRINT  939.I.P(2*J)#RH</.J)»E(2»J) 

960  continue 

liltl 

PUNCH  910.1, II.PTUL 

print  «ii.i*ii.ptul 

9  6 1  continue 
PRINT  933 
60  CONTINUE 
STUP 
end 


subroutine  dumhcn) 

common /a a  1 n/u( 2. 1000 ) .*(2. 1000) »u(2»looo) »C(  ?» looo> ,rh( 2. 1 0 

100)  »£(  2, 1000) »P(  ?» 1 0 0 0 ) 

CQMHON/HEFL/TREF »t 

COMMUN/T IMUU/DT.UU1  ,UU^»  I  ,XMU 
C0MMUN/N0C0N/IS1.IS2.  IS>3.  IS«.  INT1.INT2.IMAX 
COMMON/NUI M/TT. TTMAX. X*Z»TMAX 
C0MM0N/GAM/GAMM(2) 

C0MMUN/C0N/TCUN.DT2 
COMMUN/NNEww/KSHOCK.KTELG.  IDT 
REWlNU  N 

WRITEC N)T.DT.TCUN 
WRITEC N)GAMMC 1 ) »GAMM(  2 ) 

WR I TE ( N ) UU  1  ,  UU2 
WHITECNJIS1.1S2.1S3. ISM, IMAX 
WRITE CN) INTI.  INT2.KTEL9.KSH0CK 
UO  9500  J* 1 » I  MAX 
WRITECNJXC  1,J),U(1»J),C(1,J) 

WRITECNJPC 1,J)*RH(1,J)*E(  l.J) 

9 s o o  Continue 

WRITECNJ1 , p  TUL 
REWIND  N 
RETURN 
End 


subroutine  reaoucn) 

COMMON/ GAIN/0C2. 1000) »  *  {  2  » 1000) »U( 2. 1000) ,C(  ?»  1000) .RM( 2, 1 0 
100)  .E(2»1000)»P(2.1000) 

COMMON/REFL/TREF.T 
COMMON/TIM UU/DT.UU1  »UU<?»  I.  XMU 
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CoMM0N/NQC0N/I$i»IS2.1*3»IS4.INTl»lNT2»IMAX 
CQMMON/NO I M/TTiTTMAX.XXZ.TMAX 
COHMON/GAH/GAMM(2) 

COHM0N/C0N/T60N»0T2#PTUL 


".1PE15.8/ 

".1PE15.8/ 

"•1PE15.8) 


COMMON/NNERW/KSHOCK.KTELG. IDT 
9Q?  FORMAT! 1 HO»"NON*OlMENSlQNAL  INPUT  DATA"/ 

(1H  ."(T)TIME  OF  THE  INITIAL  TIME  LINE 
llH  » " C  DT  )  T I  ME  INCREMENT  TO  THE  FIRST  TIME  LINE 
(1H  ."(TMAXIMAXIMUM  RUN  TIME 
VOH  FORMAK  1H0#"GAMMA  FDR  REGION  ONE  ".1PE15.8/ 

(lH  ."GAmMA  FQH  REGION  T*0  ".1PE15.8)  clArtTv  *  1Dri<; 

909  FoRMAT(1HO."(UU1)LEFT  TRAVELING  SHOCK  VELOCITY  ".PEIS. 8/ 
(lH  UU2)RIGHT  TRAVELING  SHOCK  VELOCITY 
9 1 J  FoRmATUHI. "PROPERTIES  OF  THE  INITIAL  TIME  LINE/ 

(  lHO ."POINT  NO.  "  »7X."*"»  1«X»’*U",UX."C".14X»  P  .13X.  RH  .1 


14X#"E") 

9t*»  F  QRMAT  (  1  H  .3X.lM»3X#6(  1PE15.8)  ) 

9?5  FUHMAT(1HO."IPTUL)SHOC*  PRESSURE  JUMP  TOLERANCE 
HEW  I  NO  N 

READ  CN)T.[)T#TCUN 


".1PE1S.8) 


Tmax*tmax/tcun 

I F ( IDT.EQ.I )0T-0T2/TC0N 

HEAU  (N)GAMHtl).GAMM(?J 

R£AU  ( N  )  UU l » UU2 

READ  (N)IS1»1S2»IS3.IS<*.IM4X 

READ  (N)INT1»INT2.KTEL«.KSH0CK 

PRINT  907.T.0T.TMAX 

print  roa.gamm! i ) »gamh! 2 > 

PRINT  909.UUl.UU? 

PRINT  913 
DO  9500  J*1»IMAX 
READ  (N)X(1.J)*U(  1 » J ) » L (  l.J) 
HEAD  (N)P(1.J).RH(1.J)*E( i.jt 
Q( 1 » J)»O.Q 

print  914»j»xu»j>»u(1»j>»c(1* 


J )  »P( 1 » J ) » RH( l.J).E(l.J) 


9500  CONTINUE 

head  tNH.PTOL 
PRINT  925.PTUL 
REMIND  N 

return 

end 


SUBROUTINE  ADDPTS(L)  uuo  1 

CoMMON/GAIN/  Q(2»1000)»X(2»10QO)»U(2»1000)»C{2» 1000).RH(?»1 

1000).  L(2. 1000)»P( 2» 1000) 

CoMMON/TIMIIU/  DT.  UU1.  UU2  *  I.  XMU 

CflMHQN/NOCON/  I S 1  *  I S2 *  I  S 3 .  I  S4 .  I N T 1 »  INT2.  IMAX 

CuMMON/7n1T/PI1.UI1.RhII.EI1»CI1»0I1.PI2.U|2.RHI2.EI2.CI?.0 
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112 

COMHON/SINPT/  PI#  UI»  KhI#  UFi  XZ»  PF»  IR1»  IR2 
Common  /  otts  /  ott 

c  aoo  fohmaK"  aou  puints  tu  the  region  right  to  main  shock”) 

flo2  F ORMAT ( 1 H  #  12#  IX#  U#  3X#  6(£15#8»  2X)>"  ADDED  PTS") 
ISPEC*IMAX-1 
IS4P1  ■  IS4  ♦  1 
IMAX  •  154  ♦  10 

IF(ABS(0T)  .LT.  0#000l  )  GO  TO  5 
c  If  (  I  *E0.  1  )  OX  ■  0 T ^ C 1 2 

I F (  I , Eo • 1 )  DX»0T/C(2#lsPEC)*50. 
1F(I.UE«2)0X"UTT/C(2»I$PEC)*50* 

C  I F(  1  .GE#  ?)  DX  *  D  T  T  ^  C I  2 
S  DX-2.0 

DO  10  J  ■  IS4P1#  IMAX 
X ( 2 #  J)  »  X(2#  J-l)  ♦  Dx 
P(2#J)*P(2»ISPEC) 

U(2# J)«U( 2# 1SPEC) 

E ( 2  #  J ) »E ( 2  » iSPLC) 

C(2» J)»C(2# 1SPLC) 

RH(2# J)*RH( 2#  ISPEC) 

0(2# J)*Q(2f ISPEC) 

10  continue 
keturn 
End 


function  ADT  (M#N) 

common /gain/  u( 2# ioooj »x( 2# loooj  .u(2»  iooo)#C(2»  iooo.rhi  ?*  1 
1000).  E(2# 1000)*P(?» 1000) 

COMMON  /TIMUU/  DT.UUl #uU2#  I »  xmu 
loi  format  < i m  ,"dt  is  negative". 2(5x»I3>> 

RA  ■  0  «  99 

DX  ■  XC 1 #M)  -  XC  1  #N) 

AT  ■  DX/(C(1»M)  ♦  C(1»N))  *  ?  *  0  *  R  A 

IF  C  AT  .IE.  0.)  PRINT  101#M.N 

A[)T  ■  AT 

RETURN 

END 


SUURUUTINE  asigpt 

C  TO  ASIGN  PTS  FOR  THE  INITIAL  LINE#  FOR  MSR  *  1  CASE 

CQMMON/G A  IN/  U ( 2 . 1 000 ) *  X < 2 #  1 000 ) » U ( 2 # I000).C(2#  1000) #HH(  2# 1 
1000)#  £(2#100Q)»P(2'1000) 

COMMON/RAWAV/  XRI 30) »UH(  30) #CR(  30 ) # RHK C  30 ) »ER(  30 ) # PR<  30 ) . HH 
1P( 30) 

COMMON  /C1AND2/  PF1»UF1»RHF1.EF1#CF1»QF1#PF2»UF2#RHF2»EF?»C 


50 


r>  o  o  o 


1F2#QF2»  XF.XS 

COMMON  /  I N I T /  Pll#UIl'RHIl»EIl#CIl»Qll»PI?#Ul2#RHT?#EI?#CI 

12.012 

COMMON/nOCON/  ISl#  IS?»  I S  3 »  IS**.  INTI#  INTP.  IMAX 
COMMON/SINPT/  PI#  UI.  *Hl#  UF#  X7»  PF#  IR1#  IK2 
COMMON/TIMUU/  DT#  UU1#  UU2#  I#  XMU 
COMMUN/SHKI/  EP 

C  I H 1  a  FIRST  PT  IN  RAw*vEl  IR2  ■  LAST  PI  IN  RAwAvE 

moo  format  uh  #ix#"i"»"  j".iox#"x"»i7x#"u*»i6x#,*c"»16x#"P"#i 

1 6  X  #  nRHn  .  16X#"E"» 10X*"L") 

no l  format  c  i h  #I2#ix»I4»3X#6(E15.8#2X)»I2) 
ho2  format  (  i h  ,¥x#"q  *  "»eis.8) 

b i o  Format ( *  initial  data  being  read  in  and  arranged  in  ordfr" 
i»  /#  "  second  shock  also  been  inserteo  at  isi  and  i 

2S2") 

B  l  1  FORMAT!**  I  s  1  a"#  1 4 «  "  IS2  *"#  14#  "  I S  3  *" »  l4»  "  I.S4 
1«"#  1 4  #  /#  "  INTI  «"»  1 4 »  "  I NT2  *"»  14#  »  IMAX  a"#  14) 

b 1 2  format  uho#"  eno  of  initial  data  arrangement  for  blast  wav 

It#  READY  FDR  STARTING  CALCULATION") 

C  CALCULATE  AVERAGE  OT  In  THE  RAREFACTION  FAN 

0  T  1  *  0. 
lSHOK  *  1 
I  R  1  R  a  IR1 
IR2P  a  IR2 
1 K 1 2  *  IR2  -  IR1 
DO  b  N»IR1#IR2 
b  X( 1 .  N)  a  XR(N) 

DO  10  Mm2  » I R2 
NmM-1 

DTI*  AOT  (  M.N)  ♦  OT 1 
tO  CONTINUE 

X  I R2  *  1 H 2  -  1 R  1 
AVDT  *  OTl/XlR? 

c  CALCULATE  OX  AND  VARIABLES  BETWEEN  CENTER  LINE  AND 

c  rarefaction  wa 

N 1  *  X( 1#  I  R 1 ) / ( C 1 1  *  A VUT  > 

C  PRINT  B16.lRl.IH2»DTI.AvnT»CIl»X(  1#IR1)#N1 

C  B 1  6  FoRMAT(////#"IRla"#I3."lR2a",I3»"QTl*"»El3.6.nAVDT*,,#£l3.6# 
C  "  C  1 1  *  "  # 

C  i  Eli. 6#"  X(l»IRl)*".El3.6."Nla",l4) 

1R1  *  N 1 ♦  1 
I  H  2  ■  I  R  1  ♦  I R  1 2 
00  15  NmIRl.IR? 

X( 1 »  N)  *  XR(  IR1P  ♦  N  •  IR1) 

P( 1#  N)  m  PR(  IR1P  ♦  N  -  IR1) 

U( 1 #  N)  *  UR (  IRlPfN-lRl  ) 

PRINT  817.N.XC 1 »N)#P(1 »N) »U( 1 »N) 

B 1 7  FORMAT!  1 H  . "N*" , I  4 » " X  (  1 # N ) a" , E 1 3 . 6 » »P< 1 # N ) a" , £  1  3 . 6 # "U(  1 » 

N)**#e1 3.6 
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C  *) 

HHU»  N)  •  HHH(  IR1P*H-IH1  ) 

E  (  1  »  N)  *  ER(  IRlP*N-lRl) 

C ( 1 »  N )  •  CRC  IHlP+N-lRl) 

15  Q( 1 ♦  N)  «  0. 

XNl  ■  N 1 

OX  ■  X(  1  ,  IR1 )  /  X  N 1 
X(  1  » 1 )«0. 

IR1M1  ■  IR1  -  1 
00  20  N«2, IR1M1 
X(  1  »N)*X(  1 »N-1  )  ♦  DX 
20  CONTINUE 

U[]  25  Na  1  #  I  R  1  M  1 
P( 1#N>  a  PI  1 
U( 1 ,  N)  a  UI1 
HHC  1 ,N)  *  RHl  1 
L (  1  *  N )  a  £11 
C(  1 ,N)  a  C  I  1 
U( 1 #N)  a  Oil 

25  CONTINUE 

C  “  INSERT  SHOCK  AT  PT  1*2. 

C  EP  «  0. 

c  Call  shukin  (  l*  2.  3.  IR2  ) 

C  (iu  TO  27 

C  LP  *  C  P (  1  *  I R  2*  1  )  -  P< 1 • IR2) )/P(l » IR2)  *  0.1 

c  call  shokini  i  ,  l .  l.  i*?> 

?7  CONTINUE 

C  COMPUTE  OX  AND  VARIABLES  BETwEEN  RAwAVE  TAIL  TO  INTFPT 

N2  *  C  X E  -  xR( lR2P))/(CEl*AVnf) 

XN2  «  N? 

IF  (N2  .LE.  1  )  (jO  TO  J? 

UX  ■  (  XE  -  XK(  I  H  2  P  )  )/%•'*? 

C  I N  T l  e  IS2  ♦  N? 

INT1*IR2*N? 

I N T 2  ■  INTI  ♦  1 

iNTlMl  .  INTI  -  1 

C  IS2P1  *  ISP  ♦  1 

IR2P1  «  IK?  ♦  1 

KINO  c  i 

C  00  30  N  a  I S2R 1  *  INTIMi 

N  «  1R2P1 

2»  IEUSHUK  .NE«  KING)  GO  TO  29 
NMl  a  N  -  1 

C  -  INSERT  SHOCK  AT  POIN^  NMl, 

C  EP  «  0. 

C  CALL  SH0KIN(  1 .  2*  3.  NMl  ) 

EP  *  0,001 

29  continue 

X( 1 »  N )  a  X(  t  »  N*  1  )  ♦  OX 
P(l»  N)  S  PEI 
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c 


ud»  n)  ■  ufi 

HH(1»  N)  •  RHF1 


E(l* 

Cd» 

ud» 

KING 

1F(N 

N  ■ 


id  •  EH 
N )  ■  CFl 
N)  ■  OF  1 
■  KING  ♦  1 
.EG*  INTIMI) 
N  ♦  1 


GO  TU  30 


GO  TO  28 

30  CONTINUE 
GO  TO  35 

32  INT 1  ■  IS2  ♦  1 

1  NT 2  *  INTI  ♦  1  _ 

compute  ox  and  variables  between 

35  CONTINUE 

X(WINTl)  «  AF 
P(UINTI)  •  PEI 
UOdNTl)  ■  UF1 
HHC 1  * INT  t  )  * 

E( I » INTI  )  ■  EE1 
C(UlNTl)  *  CF1 
Q( 1* INTI  )  a  0F1 


Xd»lNT2) 
Pd*INT2)« 
U(  1*  INT2) 
HHll* INT2) 
E(1»INT2  ) 
C(1#INT2) 
0(1*1  NT2) 
N3  *  (XS  " 


■  XF 
PF  2 
«  UF2 
•  HHF2 
■  EF2 
«  CF2 
*  UF2 

XF)  /<CF2*AVUT> 


XN3  «  N3 

IF  (NJ.LE.  1)  GO  TO  4? 
OX  •  cxs  *  xF)/XN3 


1  S3  *  INT2  ♦  N3 

ISA  *  IS3  ♦  1 

INT2P1  a  INT2  ♦  l 

IS3MI  ■  I S 3  •  1 

OQ  40  N*INT2P1 » IS3M1 

X( 1 » N  )  •  X(1#N-1)  ♦  DA 

P( 1  *  N  )  a  PF2 

U( 1 *N)  a  UF2 

HH( 1 »  N )  a  RhF2 

E( 1 »N)  a  EF? 

C( 1»N)  a  CF? 

Q (  1  *  N  )  a  OF? 

40  CONTINUE 
GO  TO  45 

42  I  S3  •  INTI  ♦  1 
IS4  «  IS3  ♦  1 

45  CONTINUE 


INTFPT  AND  shqckfpt 
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X( 1 » I  S3  )  a  xS 
P(l»l53)  ■  PF2 
U(  U  I  S3  )  ■  UF2 
HH( 1  *  I  S3 )  *  HHF2 
E( 1,  IS3)a  EF2 
C  (  1 ,  I  S3 )  *  CF2 
U  (  1  *  I  S3 )  >  U  F  2 
X(1»1S4)  ■  XS 
MlUIS4)  ■  P  1 2 
0(1, IS4)  a  U 1 2 
HH  (  1  ,  I  S  4  )  a  RH  I  2 
E ( 1 » I S4  )  aEI2 
C  (  1  » I S4 )  a  C12 
U(  1*  IS4)  a  u!2 

C  ADD  POINTS  IN  REGION  2  (  1 0  PTS) 

IS4P1  a  IS4  ♦  1 
IMAX  a  IS4  ♦  10 
DX  a  C 1 2  *  AVL)T 
00  50  N  a  IS4PI » IHAX 
X( 1 ,N)  a  X( 1 »N-1  )  ♦  OX 

P(  l.N)  a  PI? 

U(  l.N)  a  UI? 

RH(1»N)  a  RH I  2 
E(l.N)  a  El? 

C(1»N)  a  Cl? 

U(  UN)  a  01? 

SO  continue 

UT  a  AVDT 

L  a  1 
1  »  1 
Ht'TURN 
EnO 


subroutine  calcut 

common /gain/  Q( 2. 1000) *x(2.iooo>  ,ut  2. 1000) ,c<  2# 1000 .phi  ?.  1 

1000).  E( ?» 10U0 ) #P( ?♦ 1000) 

COMMON  /NOCON/  IS1»IS?MS3.IS4. INTI. INT2. IMAX 
COMMQN/TIMUU/  OT.  UU1.  UU2 .  I.XMU 
C  100  FORMAT  ( 1 H  »  1  OX  » "OT  a  ",  E15.8) 

SMDT  a  100. 

I F  (  I  NT  1  ,EQ.  0)  GO  TO 
1)0  10  K  a  ?,1S1 
XAOT  a  AOT(K.K-l) 

10  IF  (XAOT  .IT.  SHOT)  SMUJ  a  XAOT 
I  S  2  P 1  a  IS?  ♦  1 
00  15  K  a  IS2P1.INT1 
X  AOl  a  AOT  (K.K-1) 

1 5  IF  (XAOT  .LT.  SMOT)  SMUT  a  XAOT 


I NT2P 1  •  I NT2  ♦  1 
00  20  K  »  InT2P1#IS3 
XADT  ■  AOT  ( K  »  K* 1 ) 

?0  IF  (XADT  .LT.  SHUT)  SMUT  ■  XADT 
IS4P1  •  IS4  ♦  1 
DO  25  K  ■  I S4P 1 . IMAX 
A AOT  ■  AOT(K.K-l) 

?5  IF  ( X AOT  .LT.  SHOT)  SMUT  •  XAUT 
00  TO  28 

26  00  2 7  K.2  *  I  MAX 
XADT*  AOT(K.K-l) 

?7  IF (  XADT  .LT.  SMDT)  SMUT*XADT 

continue 

OT  ■  SMUT 

C  PRINT  100. DT 
return 

END 


subroutine  centptil.  id  12.  1 3 > 

CqMMON/GAIN/  Q( 2. 1000) »X(2» 1000) .U(2. 1000)  »C( 2.1000) »NH(  2. 1 
1000).  E<  2. 1000)  »P<2»1000) 

COMMON/TIMUU/  or.  UU1.  UU2.  I.  XMU 

100  FoRMATCIH  .  12.  IX.  1 4 »  3X»  6(E15.B>  2X).  12."  CENTER") 

101  F 0RMAT ( 1 1 X  .  "  NOT  CONVERGE") 

c  102  Format ( ih  ."  check  print  k*  "»m 

C20()0  FuRMAT(6X.  «  X  A" »  10X.  "UA"»  10X.  "PA".  10X.  "RHA".  9X. 

C  "EA",/, 

C  1  1H  ,  5EU.M) 

XINPtVl.  UV.  OX.  OY)  •  VI  ♦  DV*DY/DX 
LIM  ■  30 
TQL  ■  0.0005 
T0L1  *  l.E-10 
K  «  1 

C  -  DEFINE  VARIABLES. 

U3  ■  0. 

X3  ■  0. 

U( 1 »  II)  *  o . 

X  (2.11)  •  o.o 
U A  ■  U( ? .  I?) 

U5  ■  U( 2 .  13) 

0X2  ■  X ( 2 #  12)  -  X( 2.  U) 

0x3  »  X(  2 »  13)  -  X ( 2 .  il) 

0X3  «  ( U4*0X3**2  -  U5*0x2**2)/(UX2*UX3*(DX3  -  0X2)) 

Pi  •  PCI.  II) 

HHl  *  RH( 1 »  II) 

Cl  ■  C(l.  II) 
u  1  «  0. 

El  *  E(  1 .  ID 
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til  ■  ti (  1 »  ID 
P2  •  P(l#  I?) 

KH2  -  RH C 1 •  12) 

C 2  ■  C(l#  1 2 ) 

U2  *  U( 1  *  I?) 

03  ■  0( 1 »  12) 

0X1  ■  X(t.  12)  -  X C 1  *  il) 

OUl  *  U( 1,  12)  -  U<  1  #  II) 

URH1  -  RH(1#  12)  -  RH  (  i  #  II) 

□  Cl  ■  C(  1*  12)  -  C( 1.  il) 

on  ■  u l #  i2>  -  E(  i  #  ii ) 

UP  1  -  p ( 1 .  12)  -  P(  1  .  Il) 

001  ■  ti( 1 •  12)  -  ti(  1*  il) 

C  “  ASSUME  PROPERTIES  AT  POINTS  A 
PA  »  (Pi  ♦  P2)/2. 

KHA  *  (RHl  ♦  KH2)/2. 

CA  *  (Cl  ♦  C2)/2. 

UA  ■  U2 
E  3  «  El 
C 3  *  C2 
RH3  ■  RHl 
P3  ■  PI 

C  “  BEGINNING  OF  ITERATION. 

10  CONTINUE 

UMC3A  *  ( U 3  ♦  UA  -  C3  "  CA)/2. 

XA  *  *UMC  3 A  *  0  T 

UA  *  X  I NP ( U 1 #  OUl .  0X1  »  XA) 

CA  ■  X 1 NP ( C 1 »  UC1.  0X1*  XA) 

Ea  «  XlNP(El»  0E1.  0X1*  XA) 

Pa  *  XlNP(Pl#  UP1#  0X1*  XA) 

RHA  *  F.0STR0(L»  EA#  PA) 

UA  ■  XlNP(01 »  OUl »  0X1  »  XA) 

KHCjA  ■  ( RH3*C  3  ♦  RHA*Ca)/2, 

UCX3A  *  (UA*CA/XA  ♦  UX3*C3)/?. 

PE  3  *  EqSTPEIL.  E  3  *  HH3) 

PEA  *  EoSTPEtL.  EA#  PA) 

PRCU3A  .  (PE3*U3/(RM3*'-3)  ♦  PE  A  *  0  A  /  (  RH  A  *  C  A  )  ) /2  . 

PR  1 3  •  (P1/(RH1**2)  ♦  P3/(RH3**2))/2. 

U 1 3  *  (Ql  ♦  U  3 ) / 2  . 

P3P*  PA-UA*RHC3A*(-( XM^-l.)  *UC X 3 A  +  PRC 0 3 A  )  *RHC3A*0T 
E3P*  ED  PR13*(RH3-RH1  )♦  013*UT 
RH3P*  E0STR0(  L#E3P#P3P) 

C3P«  EOSTCO  (L#E3P»RM3P.P3P) 

IF  (ABS((P3P-P3)/P3)-  TnL)?l#2l»A0 
?1  IF  (ABS(  CE3P-E3)/E3>-  f0L)22»22»«0 
?2  IF  (ABS(  (RH3P-RH3)/RH3)  -  T0L)23#23#4O 
23  IF  (AdS(  (C3P-C3J/C3)-  1 Ql)  50»50#40 

ao  continue 

K«  K+l 

P3«  ( P3  +  P3P )/2  » 
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non 


E3«  (E3tE3P)/2. 

RH3«  CRH34Rh3P)/2. 

C  3*  (C3*C3P>/2. 

PRINT  2000*  X  A  #  UA»  Pa#  RHA#  EA 

PRINT  100»I«I1*X(1*I1)#U3»C3.P3»RH3»E3»L 

print  io?#k 

IFIK.LT.LIM)  GO  TO  45 
PRINT  101 
GO  TO  50 
45  CONTINUE 
GQ  TO  10 

so  Continue 

X(2#  11)  ■  0. 

U(  2#  11)  •  0. 

C  ( 2#  ID  *  C3P 
Hh(2#  11)  ■  RH3P 
£( 2#  11)  «  E3P 
P( 2 »  11)  *  P3P 
RETURN 
eno 


functiun  EQSTCQCl#  E*  «R#  P) 
C0MM0N/GAM/GAMMC2) 

101  FoHMATdHO#"  RH*"  *  E  1 5  •  8  # "  P*"#E15.B> 

I F ( RH  .LT.O.)  PRINT  lOl.RR.P 
IF(H  .LT.  0.)  PRINT  iOl.RH.P 
GO  TO  ( 1 0 »  20  )  #  L 
10  CONTINUE 

C  FUNCTION  TO  calculate  speed  of  sound 

GAHMA«GAMR( 1 ) 

C2*GAMHA*P/RH 
if  <  C2  .LE.  0  )  C?«  “C2 
EOSTCO  a  SQRTCC2) 

RETURN 
20  CONTINUE 

GAMHA*GAMM( ?) 

C2«(iAMMA*P/RH 
IF  {  C2  .LE.  0  )  C2*  *C2 
EQSTC0«S0RT( C2) 
return 
End 


function  eqsTeoil.  rh.  p) 
COHHON/r,AH/GAMM(  2) 

GO  TO  ( 10.20)  #L 
10  CONTINUE 
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GAMHAaGAHM( 1 ) 
tuSTtfl  a  P/( (GAMMA  - 

return 
20  CONTINUE 

GAMMA*CAHMC  ?) 
toSTEU*P/((GAHHA-l  ,0)*RH) 
return 
CnO 


function  eqstho  (l»e*p> 

C0MMUN/GAM/(iAHM(2) 

GO  TO  (10(2O)«L 
10  CONTlNUf 

FUNCTION  To  CALCULATE  DENSITY  FROM  Eu.  STATE 
U  AMMA*(jAMM( 1 ) 

EoSTRo«P/((GAMHm-1.U)*L) 

return 
20  continue 

GAMMA«GAMM( ?) 
toSTHO*P/(  (GAMMA-l.O)**-) 

return 

END 


►UNCTION  EOSTPE(L.E#RH> 

COMMON/GAM/GAMM(  i>) 

GO  To  ( 1 0  *  20 ) *  L 

FUNCTION  to  CAL.  DERIVATIVE  of  p  m.R.  to  e  for  RH  CONSUnI 

10  continue 

GAHMAaiiAHMt  I ) 

EOGTPE«(GAMMA-1.0)*RH 

return 
?o  continue 

GAHMAbGAMM( ?) 
taSTPE«(GAMMA-l .0)*RH 
return 
end 


functiun  EOSTPOIL.  E*  *H) 

COMMON/GAH/GAMH( 2) 

FUNCTION  To  CALCULATE  PRESSURE  FROM  EM.  STATE 
GQ  TO  ( 1 0  *  20 ) *  L 
10  CONTINUE 

GaHHA*GAMM( ] ) 

Eqstpu  ■  (gamma  -  i.o)*rh#f 
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IF  (EflSTPQ  *  t T »  0,  )  EttSTPQ  ■  0. 
RETURN 

■?0  continue 

GAMMA*GAMM( ? ) 

EaSTPO  -(GAMMA  •  1 .0)*RH*E 
IFIEUSTPO  ,LT.  0.)  EQ*TPQ*  0. 

return 

end 


function  eqstpr  (l.e»r*) 

CUMMON/GAM/GAMM{ 2) 

GO  TO  (10*20)»L 
10  CONTINUE 

C  FUNCTION  To  CAL.  DER.  OF  P  tf.R.  TO  RH  FOR  E  CONSTANT 

GAMMA*GAMM(  1  ) 

EflSTPR  *( GAMMA" l «  0  )  *E 
RETURN 

?o  continue 

GAMMA*GaHM(?) 

EqSTPR«(GAMMA-1 .0)*E 

return 

EnO 


SUBROUTINE  CALENGCXMA*) 

common/gain/  Q<2*iooo)*x(2*ioo0).u<2#iQoo>»C(2»iooO)»RH(?»i 

1000).  E(2»1000)»P(p»1000) 

dimension  SUM1(  1000) »sOm2(  1000)#SUM3(  1000) 

DIMENSION  FF(?.500)  »G'*(2.'500)»MH(2#500) 

101  FqRMA  T  (  1H  »"  KEN  EN£R^Y****E15i8»**  INT  ENERGY****  E15.8. 

I  "  TOTOL  ENERGY****  El5«8»  "IMAX*"*  15) 

102  FORMaT( 1 M  *  3E 15 • 8  ) 

1000  F0RMAT(4E15.8) 

UO*0.77fllE  OA 
U0*0 . 807263E  OA 
I*  2 
1  ■  1 

C  *****  0*7781  M/SEC.  ****  RHC J*2 * 237  GRAM/(CM)**3 

C  ********* 

D*7,7aiE  03 

C  FOR  0  IN  M/SEC. TOTOL  ENEH^Y  IN  C AL . /GRAM » THEN  CkEN*0. 239E-03* 

C  C1NT*1 .0 

CkEN*0.P39E-03 
C INI « 1 -EO 
U02*  UO  **? 

XX2*X( I  *  1  )**2 
UU2*  U ( 1 , 1 )**2 


FF(  1,1)*  E  ( 1 » 1 ) *KH(  I*  1)*XX2 

(iG(  I  » 1  )■  0,5£0*U02*RHl  I # 1 ) *UU2*XX2 

HH( 1 0 1  )»  RHl I » 1 )*  XX? 

SUM  1 1 1 )*0.E0 
SUM21 1 )«0.E0 
SUM3(  1  )  * 0  ,  FO 
J*2 

10  JM 1 ■ J“ 1 
1  1  Ox*X( I » J )"X( i » JH1 ) 

IF(DX  ,UT,  0.1E-04)  on  TO  12 
J*  J*  1 
GO  TU  11 

1?  XX2*  X( I » J)**? 

UU2*u ( I »  J ) «*2 

FF (  I »  J  )  ■  E(l.J)  •  Rh(1,J)  *  XX2 
GG (  I  *  J  )  ■  O.SEO*UO?*RH^ 1  * J) «UU2*XX? 

HH(  I  »  J  )  ■  RHl I  * J)  •  X  X  ^ 

suMiiJi*  suni(jMi)  ♦  o.seo*(FK  i »  jmi  >  ♦  ffu»j))  *u* 

SUM21J)*  SUH?(JM1)  «  U,«>E0*( GGl I  * JM1 )  ♦  GG(I»J>)  *DX 

SUM3CJ|»  SUHJ(JMI)  ♦  0,SE0*(  HHC  I  •  JMl  )  ♦  MH(I*J)>  *l)< 

PRINT  102* SUM l l J) »SUM2(  J) ,  SUM31 J) 

I F (  AhSCXf  I.J)“XMAx)  ,LT.  0  •  1 E “04  )  GO  TO  ?0 
J*  J  +  1 
GO  TO  10 

?0  SUMlNT«SUMl(j)/bUM3(j)  *C I  NT 
SUMKENaSUM?U)/SUM3(  J>  *CKtN 
I M AX* J 

ENERGY*  SUM  I N  T ♦  SUMKFN 
RETURN 
END 

SUBROUTINE  GNPSHA  (  l»  II*  12*  I3»  X4  ) 

COrtMUO/GAlN/  U(?»1000) »X(2» 1000) *U(2* 1000) *C( 2* 1000) *RH( ? . 1 
1000).  E(?» 1000) *P< ?» 1000) 

COMMON  /SHK4A  /  U 4  A  * C 4 * , RH 4  A  * E 4 A * P 4 A  *  X 4 A  *  0 4 A 
COMMON/T1MIJU/  OT  *  UU1.  UU2  *  I*  XMU 

U I  NR ( V2  »  DV1»  0V2*  Y)  *  V2  ♦  IDV2*DX1**2  ♦  OVl  *0X2*  *?  )  *  Y/ (  It 
1X1*0X2*  (0X1  ♦  OX?))  ♦  («DV1*UX2  ♦  0  V2*l>X  1  )  *  Y*  *  2  /  (  0  X  1  *  0  X  ? 

2  * ( OXl  ♦  0X2)) 

luoo  Format  <  i h  , "point  4A  ores  not  converge”) 
c i oo i  format  ( ih  »S(Eii.4)) 

cioo2  format  uho."  P4P  uap  F4P  rh4P  cap”) 
c i o o 3  format ( i h  ,"xa  *  "*eis»b."  xb  ■  "»Ei5,8."xc  *  "*ei,>.h."ua 

c  *  "* 

C  *  E  1  5  •  8  ) 

C 2 0 0 0  F QRM A  T ( "  A  "#  6(E1S.8»  2X)) 

C 2 0 0 1  FORMATC"  H  "#  6(E15«B»  2X)) 

C2U02  F  0  R  M  A  T  ( **  C  "*  fo(ElS.8'  ?X)) 

C  PRINT  1002 

T0L  «  0,0008 
ToL 1*1 , F “20 


K  ■  1 

C  -  INITIAL  ASSUMPTIONS  FOR  PT  4. 

U4  ■  UCl,  I?) 

CH  m  CC1,  I?) 

RH4  ■  RH  (  1.  12) 

£4  ■  £ ( 1  *  I?) 

P 4  ■  P(l.  I?) 

UA  ■  U(l*  13) 

CA  ■  C(  1#  I  3 ) 

UP  ■  U(  1 .  ID 
Cfl  ■  CCD  ID 
UC  -  U(l.  I?) 

cc  *  CCl.  I?) 

C  -  DEFINE  VAN  I  ARLES • 

NH2  •  RHCD  12) 

C?  *  C( 1 i  I?) 

U?  »  U (  1  *  1?  ) 

P2  *  PCD  12) 

£2  ■  £( D  I?) 

02  *  0(1.  12) 

Oxl  «  X  (  D  T  2  )  -  XC  D  ID 

DUl  *  U( D  12)  -  U( 1»  1 1  ) 

UC1  «  C(  D  12)  -  CCD  U) 

0P1*P( 1 • I  2 ) -PC  D I  1 ) 

UHH1  *  HMC  D  12)  -  RHC  I*  ID 
UE1  *  ECD  12)  “  EC  D  M) 

Dfll  ■  U  (  D  T  2 )  •  OCD  Al) 

0X2  *  X ( D  13)  -  XC  D  A?) 

0U2  ■  U( 1,  13}  -  U(  D  A 2 ) 

DC 2  ■  C(l»  13)  •  CCD  A?) 

0P2«PC  l.I3)-PC  DI2) 

UNH2  a  RHC  D  13)  -  HHC  A  .  12) 

UE2  ■  E(  D  13)  •  ECD  A?) 

002  ■  OCD  13)  -  OCD  A?) 

C  “  ES1IMATE  POSITIONS  FUr  PTS  A  »  B*  AND  C» 

io  continue 

XA  ■  X4  -  (1)4  •  C 4  ♦  U*  •  CA)/2.*DT 

XB  ■  X4  -  (U4  ♦  C4  ♦  UO  ♦  CB  )  / 2 »  *DT 

XC  ■  X4  -  (U4  ♦  UC)/2«*nT 
OXA  *  XA  -  XC  D  12) 

UXB  ■  XH  -  XC  D  12) 

0 X C ■  XC  -  X  C  1  .12  ) 

C  PRINT  1003  .XA.XH.XC.U* 

PA*OInP{P2#I1PDOP2»DXA) 

£  A  ■  OlNPCE?#  D£D  DE?»  OXA) 
RHA«EoSTRQCL»EA»PA) 

C  PRINT  2000.  PA.  EAt  RH*  »  CA.  UA.  XA 

IF  CRHA  .LT.  0,  )  RHA«-0.5D0*HHA 
CA  ■  EQSTCQCL#EA»RHA.P*) 

UA  ■  OlNPCU?#  DUl.  DU2»  DXA  ) 
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c 


c 


c 


PEA  ■  EaSTPEU*  E  A  #  RHA) 

oa  ■  uInpiq?*  dui#  no?*  nxA) 

P8*glNP(  P?#DP1#0P2.0XU) 

EB  ■  U  1NP( E?  *  DEI#  DE?»  OXB) 

RmB*EuSTRQ( L#EB*PB) 

PRINT  2001  *  Ptf *  EB#  RHB*  CB#  UB#  XB 
IF  (RHB  .LT.  0.  )  RHB«"0.5D0*RHB 
CH  ■  EQStCO  (L.EB.RHB.Pfl) 

UH  a  UlNPCU?*  OUI#  DU?  *  OXB) 

PEB  a  EoSTPEiU#  EB#  RHU) 

UB  ■  UlNPlQ?*  OUI#  0Q?»  OXB) 

EC  *  UlNP(E?#  Oil#  D£2»  OXC) 

PC  «  UlNP  (  P2»UPl»l)P2»nXC) 

RHC  a  EUSTRO  (L#EC#PC) 

UC  «  ulNPCO?#  DU  1  *  002  *  OXC) 

P£4  a  EQSTPE  (  L#  E4  #  RH4  ) 

04  *  UC 

-  CALCULATE  CUEFF.  IN  (  HE  CHARAC.  EO.. 

KHC4A  O  (  R H 4 * C <1  ♦  RHA*Ua)/?. 

HHC4B  a  (HH4*C4  ♦  RHH*t“B)/?« 

UCX4B  a  (U4*C4/X4  ♦  UH*CB/XB)/2# 

UCX4A  a  (U4*C4/X4  ♦  U  A  *  C  A  /  X  A  )  /  2  • 

P  RC04B  a  ( Pt4*U4/ ( RH4*C4 )  ♦  PEB*QB/(HHB*CB) )/?• 
P  RCQ4A  a  ( pE 4*U4 / ( RH4 * C 4 )  ♦  P£ A *0 A / ( KHA*C A ) ) / 2 « 
PEQ4A  a  ( PE  4  *U4  ♦  PEA*Ua)/2# 

P  E  0  4  U  ■  ( PE  4  *  U  4  ♦  PEBa^B)/?* 

PRH4C  ■  (P4/RH4**2  ♦  PU /RHC * * 2 ) / 2 • 


U  C  4  a  Qc 

P4P  a  (PH/RHC4H  ♦  PA/RPC4A 
I  4  UC  X4  A )  ♦  P  RUQ4B 


UB  -  UA  ♦  C(XMU  -  ?,.)*(UC*UH 
P  RC04A)*DT)/( 1  * / R H C 4 A  ♦  1*/Rh 


2C4B) 

U4P  a  (PH 
1RHC4B  ♦ 


.  pa  ♦  RHC4B*U8  ♦  RHC4A*UA  ♦  (-(XMU  ■  1.)*  UCX4ri* 
PEU4H  ♦  (XMU  -  1 . )*UCX4A*RHC4A  -  PEU4 A ) * 0 T ) / ( H *C 4 


2b  ♦  RHC4A) 

E4P  a  EC  ♦  PKH4C*(RH4  “  RH2)  ♦  QC4*DT 


HH4P  a  EQSTRU(L#  E4P.  ^4P) 

PRINT  1001 #P4P#U4P»E4P»RH4P#C4P 
IEIRH4P  ,lT.O.  )  RH4Pa“0.500*RH4P 

C«P  ■  EUSTCOCL#  £  4P  #  RH4P#  P4P) 

PRINT  100l#P4P*U4P#E«P'RH4P#C4P 

IE  (  AOS  (P4P).LE.  T0Ll)GU  TD  15 
IF  (  MBS  ( ( P4P*P4 ) /P4P ) • GT «  TUL)  GO  TU  20 
1  *3  IF  (  AHS  (U4P)  #LE  ♦  T0L 1  )  GO  TO  IB 

it  (  ABS  (<U4P  •  U4)/U4P).GT.  TUL  )  GU  TO  20 
Id  IF  ( ABS(E4P)  iLE.  TOLl)  GO  TO  20 

IE  (  ABS  ( ( E4P-E4 ) /E4P >  #GT«  TOL)  GU  TO  20 

GU  TU  30 

?0  I E ( K  #GE.  20)  GU  TO  35 
K  *  K  ♦  1 
U4a(U4p4U4)*0.5 
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o  c-> 


C4*(C4P+C4)*0,5 
P4*(P4P»P4)*Q»5 
HH4*(RH4P*RH4)*0.5 
£4«(E4P*E4)*0.5 
GO  TO  10 
ib  PRINT  1000 
30  CONTINUE 
X4A  •  X4 
UaA  ■  U4P 
C4A  «  C4P 
KH4A  ■  RH4P 
£ 4 A  -  E4P 
P4 A  *  P4P 
U4 A  *  04 


c 

PRINT 

2000.  Pa. 

E  A » 

RH*. 

CA. 

UA. 

XA 

c 

print 

2001.  pb. 

EB. 

rhu. 

CB. 

UB. 

XB 

c 

print 

2002.  PC. 

EC. 

RHG, 

CC. 

uc . 

XC 

c 

PRINT 

1003  .XA. 

XH.XC 

.  U  A 

RETURN 

ENU 


subroutine  gnpshbil.  H#  12»  13.  X4»  ms> 

COM  MON  / GAIN/  U(  2 »  1  000  )  *  X  ( 2 »  1  000  )  .  U(  2  »  1 000  )  . C  (  2  *  1 000  )  .  RH  (  ?.  »  1 
1000).  E(?»1000)»P(2*1000> 

COMMON  /  SHK4H  /  U4B . C* B » RH4B . E4B » P4B » X4B . G4B  »UUU 

CUMMON/TIMUU/  0  T  #  UU1.  UU2 »  I#  XMU 

common/gnpsbm/  XB 

loO  E  URMaT  ( ftX  »  "XB".  10X»  "UB"»10X»  "PB"  » 1  OX  »  "RHB"  »  9X  ,  "E  fl- ,  /  » 

1H  » 

1  5E11.4) 

C  1 0 1  FoRMAT(6X»  "P4B  *"»  Ell. 4#  5X»  "U4B  »"»  Ell.4) 

X I NP ( V 1 .  QV»  UX.  DY)  «  VI  ♦  DV*DY/DX 
C  M  MS»l*RlGHT  RUNNING  SMOCKl  MS  *  2.LEFT  RUNNING  SHOCK. 

SIGN  «  1  , 

I F ( MS  .EQ.  1)  SIGN  ■  -1. 

X 2  *  X(l»  I?) 

U1  ■  0(1.  I?) 

Pi  *  p(l»  I?) 
ui  ■  uu,  1?) 

El  ■  E( 1 »  I?) 

OxB  ■  XB  -  x(l*  12) 

IF ( XB  .LT.  X 2 )  GO  TO  1^ 

UX  «  X( 1 »  I  3 )  -  X(  1.  I  2  ) 

OQl  *  0(1.  H)  -  0(1.  1 2 ) 

0P1  ■  P( 1 .  Ii)  •  P( 1  *  1 2 ) 

DEI  ■  E(  1  *  I  3 )  -  £(,  1.  1?) 

DU  1  «  U(  1 .  13)  -  U(  1 .  1 2 ) 

GO  TO  20 


£L  7 


1C 


10  CONTINUE 

Ux  ■  XCl.  tl)  -  XU#  \'i) 

DPI  ■  P<  1.  ID  -  PCI.  *2) 

U£1  ■  E  C  1 »  ID  “  E(h  1?) 

OUl  -  U(  1  #  ID  -  U(l,  X  2  > 

UQl  ■  U(  1 »  ID  -  till#  T  2 ) 

?0  CONTlNUt 

PH  *  XlNP(Pl*  DPI#  OX.  OXB) 

EH  a  XInPCEI.  ntl.  OX.  OXH) 

KhB  ■  EQSTROIL.  Ed.  PH) 

CB  «  EOSTCOC  L »  td.  RHH »  PB) 

UH  a  XlNPiUl#  OUl,  OX.  OXB) 

UB  »  XlNPCQl,  OUl.  DX*  OXB) 

PEB  *  EtiSTPEtL.  EH,  RHH) 

HhC  a  RhB*CB 
UCX  *  UB*CH/XH 
PRCQ  «  PEd*Od/RHC 
U4B  «  O* 

P4d  *  PB  *CSlGN*CU4B  -  Ufl)  ♦  ( “  (  X  MU  -  l.)*UCX  ♦  PHCGDOT)  * 
1  RHC 

IF  CP4B  •  LF .  0.  )  P4B“-P48 
L  PRINT  100.  XH.  UH,  PH.  RHB .  EH 

C  PRINT  101,  P4H  »  U4R 

RETURN 
EnO 


SUBROUTINE  f,NPT(L.  II,  12.  II.  KQ  ) 

CQMMON/GAIN/  U( 2. 1000) »X( 2 .1000) ,U(  2. 1000) ,C( 2. 10U0) ,RH( ?, 1 
1000).  E(?» 1000), P(?* 1000) 

CALCULATE  PROPERTIES  AT  GENERAL  POINTS* 

COMMUN/T I MUU/  OT.  UU1.  UU2,  I.XMU 
REAL  LINP 

12.  IX,  14*  3  X  »  6 ( E 1 5 . 8  ,  2X).  12.  2X.  2HH1) 

12.  IX,  I  4  *  3X.  6(  E  15  •  8 ,  2X  )  »  12.  2X.  ?HM?) 

12,  IX.  14.  3  X ,  ft ( E 1 5  *  8  »  2  X  )  .  12) 

10X»  "NOT  DIVERGE") 

"PT  Ml  CALU,  PT  H  LIES  ON  LEFT  RUNNING  SHUCK. 


1()0  FORMATCIH 
lOl  FuRMAHIH 
1 0  3  f ORMAT(  1  H 

1  (-ormatcih 
JoO  FORMAT! 1H 
1") 

JOl  FoRMATCIH 
1 ) 

400  FoRMATUH 
4  0 1  FORMATCIH 
1") 


"PT  M2  CALU.  PT  A  LIES  ON  RIGHT  RUNNING  SHUCK" 


"XH  *  ".ED. 8,"  0T1  ■ 

"1ST  CHARAC ‘ERISTIC  DOES 


".E15.8) 

NOT  INTERSECT 


The  shuck 


402  FORMATCIH 
1") 

C  40 3  FORMATCIH 
C  5oO  F URMATC  1H 
C  50 1  F ORMATC  1H 


"2N0  CHARACTERISTIC  DOES  NOT  INTERSECT  THE  SHOCK 

"  XA  ■  ",Ei5,8»"  DT2  •  " » E  1 5 * 8 ) 

"PUINT  A  VARIABLES  FOLLOW  ") 

"PUINT  H  VARIABLES  FOLLOW  " > 
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uuuuuuuuu 


c  5o2  format  mm  »**  x  u  c  p  rh  e  «i 

510  F  llRMAT  ( "  XB  NOT  CONVERGE  IN  GNPT") 

5?o  formats*  xa  not  converge  in  gnpt**> 

510  forma t ( *•  from  gnpt  we  have  U3  eis.8#  »  ou4  ■«.  eis.b. 

1  "  0X4  ■  E15.8.  «  X-*  ■  E 1 5 « 8  »  /»  «  Cl  ■"»  E15.B.  " 

OCA  ■ 

2  E 15 « 8  ,  M  XA  ■  "»  £  1 5 • 8  ) 

*  KG  -  II  GENERAL  POINT  ADJACENT  TO  LEFT  RUNNING  SHOCK. 

-  FIRST  CHARACTERISTIC  INTERSECT  WITH  SHOCK. 

-  KG  *  21  GENERAL  POINT  ADJACENT  TO  RIGHT  RUNNING  SHOCK. 

“  SECONO  CHARACTERISTIC  INTERSECT  WITH  SHOCK. 

-  KO  •  31  REGULAR  GENERAL  POINT. 

LINP  (  vl.  OV.  OX.  OY)  ■  VI  ♦  DV  *DY  /  DX 

GINPt  V?.  OVl.  UV2.  0Y>  ■  V2  ♦  ( DV2*DX1 **2*DV1 *UX2**2 ) *  OY 
1 *  I  0Xl*DX2*(0Xl  *  0X2 )  )  ♦  (-0Vl*0X2  ♦  0V2*DXl  )*(  OY 

2  ) **2/( OX  1 *0X2  *(0X1  ♦  0X2)) 

MT  *  20 
HnIT  *  A 

T0L1  «  l.E-10 
TOL  *  0.0005 
KP  ■  KG 

GO  TO  (201.  202.  200).  KQ 
20 1  PRINT  300 
GO  TU  200 

20?  print  301 
200  CONTINUE 
C  PRINT  502 

U1  ■  UC1,  II) 

El  *  E(l.  ID 
Rh l  *  hh(  1 .  in 
oi  ■  Gd,  II) 

Cl  «  C ( 1 .  ID 
Pi  ■  PCI.  ID 
U?  lUll .  12) 

E2  «  E ( 1 ,  I?) 

RH2  «  RH ( 1 »  12) 

02  »  U( 1 .  12) 

P2  «  P(  1 .  12) 

XI  ■  X(  1 .  ID 
X?  .1  XC 1  •  12) 

X 3  •  X( 1 .  13) 

C 3  ■  C (  1  ,  13) 

U 3  ■  U  (1.13) 

E3  ■  E(  1 .  13) 

P 3  *  P( 1.  13) 

03  »  0(1,  13) 

K  ■  1 

KM  1  *  3 

UA  *  U( 1 ,  1 2 ) 

oxi  ■  xc  1 ,  i2)  -  x(i»  Id 
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c  this 


OUl 

• 

U  (  1  , 

12) 

* 

ud,  in 

uci 

m 

C  (  1 , 

12) 

m 

cd,  in 

uei 

a 

Ed, 

12) 

- 

Ed,  in 

ORH  1 

■  RH  (  1 

,  12) 

-  RH (  1, 

11) 

DPI 

a 

P(  1, 

12) 

- 

pci,  in 

DO  1 

a 

0(1, 

12) 

- 

ud,  in 

0X2 

a 

X(  1  , 

13) 

m 

X(l,  i 2 ) 

0C2«C( 1 ,13) 

•C  (  1 

, 

12) 

0U2 

a 

U(  1  • 

13) 

- 

Ud,  i  2 ) 

0E2 

a 

E(  1 , 

13) 

- 

Ed*  i?) 

IS 

the  point  of 

THE  INFAMOUS 

0P2 

a 

P(l» 

13) 

- 

P(l,  i?) 

0RH2 

«  RH(  1 

,  13) 

“  RH( 1 » 

12) 

002 

a 

0  (  1  . 

13) 

• 

0(1,  i?) 

GO  TO 

(?,  4 

,  10) 

*  KQ 

H 

LIES  ON 

The 

LEFT  SHOCK, 

L)X  3 

a 

X(  1  , 

1 1 ) 

m 

X( 2,  il) 

0U3 

a 

U(  1  , 

1 1 ) 

• 

u( 2*  in 

OE  3 

a 

E  (  1 , 

1 1 ) 

m 

E ( 2,  il) 

OP  3 

a 

P<  1 , 

1 1 ) 

- 

P ( 2  *  il) 

ORH  3 

*  RH  (  1 

.  ID 

-  RHd, 

in 

00  3 

a 

Q(  1  , 

I  1  ) 

• 

0(2,  il) 

0C3 

a 

C(  1  , 

I  1  ) 

- 

C ( 2,  il) 

X  l  P 

a 

X( 2, 11  ) 

X  4  * 

X2  ♦  U  ?  *  0  T 

UTl 

a 

OT/2. 

XH  *  (XI  ♦  XlP)/2. 

KCOUNT  a  1 

HI  ■  U1 fCl-( XI *X1P)/0X3*(DU3*DC 3DDX3/DT 
d2*XlP*CUl+Ci«Xl*(0U3*UC3)/0X3+X«/XlP*DX3/0T) 
H  3**  l./OX3*(UUJ»OC3) 

s  continue 


FP«l,-2.*d3*Xd/dl 
F  a  XH  -  (  B2  ♦  33  *  X ^  **2)/dl 


XHP*Xd"F/FP 

IF  (  AdSdXRP  -  Xd)/X«P)  ,LT  ,  TUL)  GO  TO  6 
KCOUNT  a  KCOUNT  ♦  1 
iFCKCUUNf  .C.T.  NIT)  GO  TO  7 
XH*  XdP 


C  PRINT  400»Xfi»UTl 

GO  TO  S 

7  print  5 i o 

C  PRINT  4Q0,  Xd,  OT1 


6  CUNTINJF 

UT IP  *  OT / ( X 1 P  -  XI)  *  (  X1P  -Xd) 
UTl  *  OT 1 P 

I F  C  L)  T  1  ,GT.  1  #  0  *  D  T  )  GO  TO  3 
KM  1  *  1 
GO  TO  10 
3  KM1  *  3 
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PRINT  401 
GO  TO  10 

c  -  a  lies  on  the  right  shock* 

4  0X4  *  X(  1*  13)  *  X  ( 2  *  *3) 

DU4  ■  U(  1 »  13)  -  11(2*  A  3 ) 

0E4  •  E( 1  *  T  3  )  -  E(2»  I  3 ) 

DP4  •  PCI*  133  “  P( 2*  a3) 

0RH4  *  RH (  1  #  13)  -  RH(^.  13) 

0Q4  ■  0(1*  13)  -  U(  2 »  A 3 > 

0C4  ■  C( 1 »  13)  -  C ( 2  *  I  3 ) 

X3P  ■  X ( 2 » I  3 ) 

X4  •  X2  ♦  U?*OT 
XA  ■  (X3P  ♦  X3)/2. 

UT2  ■  OT  /  2* 

C  PRINT  530*  U3*  UU4  *  0X4,  X3*  C3»  0C4*  X4 

KCOUNT  •  1 

A 1  a  U3-C3-(X3*X3P)/DX4*{DU4-0C4)fDX4/DT 

A2  *  X3P*(U3-C3-X3*(DU4«DC4)/DX4*X4/X3P*DX4/DT) 

A3**  1 ./DX4*(UU4*DC4> 

9  continue 

Fp  *  1  •  *>2  *  *A3*  XA/A 1 
F  «  XA  -  (  A2*A3*XA**2>/A1 
XAP  *  XA  -  F/FP 

IF  (  A8S((XAP  *  XA)/XAP  )  ,LT.  TOL  )  GO  TO  11 
KCOUNT  a  KCOUNT  ♦  1 
IF(KCUUNT  ,GT.  NIT)  GO  TO  12 
XA  «  XAP 

C  PRINT  403*  XA.  UT2 

GO  TO  9 
12  PRINT  520 

c  print  403*  xa»  UT2 

ll  continue 

0T2P  a  OT  /  t  X3P  -X 3 )  *  (  X3P  -XA) 

UT2  ■  DT2P 

IF ( DT2  .GT.  1.0*0T)  GO  TO  8 
KMl  a  2 
GO  TO  10 
8  KMl  a  3 
PRINT  402 

C  *  BEGINNING  uf  iteration  LOOP* 

10  CONTINUE 

X4  a  X2  ♦  U2  *  OT 
KP  a  KMl 

GO  TO  (15.  ?Q.  25).  KP 
15  CONTINUE 
DT2  a  OT 
UX8  *  XB  -  Xl 

U8  *  LlNP(Ul*  0U3*  0X3*  OXB) 

Eb  ■  LlNP(El»  OE 3*  0X3*  DXB ) 

PB  ■  LlNP(Pl*  OP 3 *  0X3»  DXB) 
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NHB  ■  EqSTRq l L  #  Ed#  PH  i 
Ch  a  LUSTCO(L.  LB  »  RHH »  P8) 

UB  *  LlNP(Ql#  DU  3  »  0X3*  OXH) 

P£B  ■  EttSTPEd.Ed.RHR) 

GO  TO  30 
?0  CONTINUE 
DTI  a  OT 
UXA  »  XA  -  X3 

UA  *  LlNP(U3*  OUR*  0X4*  OXA) 

LA  *  LINPIE3*  OL  4  *  0X4*  OXA) 

PA  *»  LlNP(P3*  DP4#  0X4 »  OXA) 

NHA  a  EttSTRQd#  E  A »  P  A  ) 

CA  ■  LQSTCttd.  LA#  HH A  *  PA) 

UA  «  L I N P ( Q 3 *  DU4*  0X4*  OXA) 

PEA  a  LOSTPE(L»EA»RHA) 

<i Cl  TU  28 
UT1  a  UT 
UT2  a  DT 
Sti  CUNllNUE 

IF  (  K  .NE.  1  )  GO  TO  ^ 6 

XH«(  X4-UT*(UUCl-Xl/0X  W(DUUOCl)  )  )/(  1 .  O^UT/OX  U  (  OU  1  fDC  1  )  ) 
(iO  TO  2y 
7b  CONTINUE 
7?  UXB  *  X H  “  X 2 

UB  *  UlN  P(U2.  UUl,  OU^.OXB) 

LB  a  ulN  P  ( F  2  *  0L1.  UE^.OXB) 

UB  a  UlN  P ( 02  *  UUl*  Dtt^.OXB) 

PR  a  UlNP(P?»  DPI*  DP2  *  0X8) 

KHB  a  Los  T Ro ( L  *  Eb#  PR> 

CH  a  LUSTCQ(L.  LU»  RHH  *  Pfl) 

PEB  a  EttSTPEd#  EH.  RHU) 

IF  (  KP  « EO •  3)  GO  TO  30 
liG  TU  40 

continue 

DXB  a  XH  -  xl 

dm  a  LInPC  III  *  0U1  »  0X1  »DXH) 

Lb  *  L 1 NP ( E 1  *  OL 1 #  0X1 »DXB) 

UB  a  LInPCQI*  OUl#  DXl'OXB) 

PR  a  LlNP( P) *  DPI #  0X1  *  0X8) 

HHB  *  CuSTRttd#  Ed#  PH) 

CR  a  EfcSTC(5(L#  LB*  RHH »  PB) 

PEB  a  EoSTPEd.  LB*  RHB) 

IF(KP  .FQ.  3)  GU  TO  30 
GO  TO  40 
30  CONTINUE 

IF  (  K  ,NL.  1)  GO  TO  3<* 

XAa(X4-DT*(U3-CJ-X3/DX^*(0U2-UC2) ))/( l*O»OT/0X2*(OU?-UC2) ) 
GO  TO  31 

32  CONTINUE 

33  OXA  a  XA  -X2 
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o  r: 


UA  «  QlN  P(U 2.  0U1.  DU^.DXA) 

Ea  »  fllH  P(E2»  DEI*  OE2.DXA) 

OA  ■  QlN  P( 02 »  DO  1 »  D&2.0XA) 

PA  ■  Q I NP( P2  »  DP  1 »  OP2*  DXA ) 

KHA  ■  EoSTROC L»  E  A »  PA> 

CA  ■  EQSTCOtk.  EA.  RHA »  PA) 

PEA  ■  EttSTPE  ( L »  EA»  RHM 
GO  TO  40 
3l  CONTINUE 

Ox A  ■  XA  -  X2 

UA  «  LIN  P  (  U2  .  0U2.  0X2. OXA  } 

Ea  *  LIN  P  (  E2 »  Q£2.  0x2.  DXA  ) 

OA  ■  LIN  P  (  02*  DO? #  Ux2»  DXA  ) 

PA  •  LlNP(P?»  DP2.  DX2*  DXA) 

HHA  ■  EqSTRQ ( L  »  EA.  PA) 

CA  »  EQSTCQtL.  EA.  RHA*  PA) 

PEA  *  EOSTPEIL.  EA.  RH#) 

4(»  continue 

IF  (  K  .NE.  1  )  GO  TO  ^5 
04  ■  02 

KHC4A  ■  RHA  *CA 

RHC4B  •  RHB  *  CB 

UCX4H  *  UB*CB/Xti 

UCX4A  «  UA*CA/XA 

P  RCQ4B  •  PE«*Qb/(RHb*^fl) 

P  RC04A  >  Pf  A*QA/  (  RHA**«  A  > 

PE04A  ■  PE A*OA 
PE04H  «  PER*«B 
GO  TO  41 

35  KHC4A  ■  (RH4*C4  ♦  RHA*0A)/2. 

HHC4H  *  (RH4*C4  ♦  RHB*^R)/2. 

UCX4B  *  (U4*C4/X4  ♦  UB*CB/XB)/2« 

UCX4A  «  (U4*C4/X4  ♦  UA*CA/XA)/2. 

P  RC04B  ■  CPEB  *OB/(RHB*CB)  ♦  PE4*Q4/(RH4*C4))/2« 

P  HCQ4A  ■  CPEA*OA/(RHA*CA)  t  PE4*Q4/ ( RH4*C4 ) )/2. 

PC04A  «  ( PEA  *0  A  ♦  PE4*04)/j», 

PE04B  ■  (PER*OB  ♦  PE4*04>/2« 

4 l  Continue 

C  PRINT  500 

PRINT  103.1 » I2.XA.UA.CA.PA.RHA.EA.L 
PRINT  501 

PRINT  10  3.1  »  l2.XB»UB.Ct),PB,RHb.EB»L 
TQLCON  ■  0.0005 

P4P  »  (PH/RHC4B  ♦  PA/R*C4A  ♦  UB  •  UA  ♦  («(XMU  •  l.)*UCX4R  ♦ 
1  P  HC04B)*DT1  ♦  ( -( X«U  -  l.)*UCX4A  ♦  P  RCU4A ) #DT2 > / ( 1 . /RH 
1 C 4B  ♦  1./RHC4A) 

U4P  ■  (PH  -  PA  ♦  RHC4B*U8  ♦  RHC4A*UA  ♦  (-(XMU  -  l.)#UCX4B*R 
1 HC  48  ♦  PEQ4B )  *0T 1  -  (-(XMU  -  1 . > *UC X4A*RHC4A  ♦  PEQ4A)*DT 
22)/(RHC4A  ♦  RHC4R) 

IF  (  K  .NE,  1  )  GO  TO  “3 
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RH4  ■  RH2 

4 5  tap  ■  E2  ♦  P2/(RH2**2)*(RH4  •  RH2 )  ♦  Q2*DT 
rh4P  *  eqstrocl.  tap.  PaP) 

IF  (  ABS  ((  RH4  -  RH4P  )/RH4)  .LT.  TOLCON)  GO  TQ  46 
C  5?B  FORMAT!  1 H  .  *  RH4a"  ♦  E  1  5  • 8 »  "  RH4Pa" » E  1 5 . 8  ) 

C  PRINT  528.RH4.WH4P 

RH4  «  (RH4P  ♦  RH4 )  /  2*0 
GO  TU  45 

43  t4P  ■  E?  ♦  (P2  ♦  P4P)/!RH4**2  ♦  RH2**2)*(RH4  -  RH2 )  ♦  (04  * 

1  02 )  /  2 »  *0T 

RH4P  ■  FOSTRU(L»  E4P .  P4P) 

46  C4P  ■  EfiSTC0(L»  E4P.  RH4P.  P4P) 

IF  (  K  .EU.  1)  GO  TO  50 

IF  (  ABS  ( ( P4P  -P4)  / P  4  p )  -  TULCON)  42.  42.  50 
42  IF  ( AdS( U4P ) »  LT  »  T0L1)  GO  TO  44 

IF  (  ABS  ((  0 4 P  -  U4)/04P)  -  TOLCON)  44.  44.  50 

44  IF  l  ABS  CCE4P  -  E 4 ) /  t«P)  -  TOLCUN)  7 0.  70.  50 
so  IF  l  K  1 G£ •  NIT  )  GO  TO  60 

IF  (  K  .  EU »  MNIT  )  GO  <0  70 
C  PRINT  103. I . I2.X4.U4P.04P.P4P.RH4P.E4P.L 

K  *  K  ♦  1 

I F ( K  .GE.  3)  oO  TU  55 

U  4  a  U4P 

04  *  C4P 

P4  *  P4P 

RH4  *  RM4P 

E 4  *  E4P 

P t  4  a  EOSTPF  (L.E4.RH4) 

GO  TU  10 

S5  U4  a  (U4  ♦  U4PJ/2. 

C  4  *  (  C  4  ♦  C  4  P  )  /  2  • 

P  4  *  IP4  ♦  P  4  P ) / 2 • 

KH4  a  ( RH4  ♦  HH4P )  /2. 

E  4  a  (  E  4  ♦  E  4  P  )  /  2  • 

PE  4  a  EQSTPE(L.  E4»  RH*M 
GO  HJ  10 

C  *  END  GF  ITERATION  LOOP. 

ho  print  i 

70  X12.  12)  a  X 4 

0(2.  12)  ■  U4P 

IF  (P4P  ,LF.  0  )  P4Pa"P4P 

IF  (RH4P  .Lt.  0)  RH4P*  -RH4P 
C(2.  12)  ■  C4P 
RHC2.  I?)  a  RH4P 
E ( 2  ♦  12)  *  E4P 

P ( 2  *  12)  *  P4P 

0(2.  12)  a  o(l,  12) 

RETURN 


SuBRUUT  J  Nt  INIOAT 
DIMENSION  A !  AO ) 

COMMflN/1  IMUU/  OT»  UU1#  UU2»  I.  XHU 

CQMHUN/lNlT/RIl#un»RHil.EIl.Cll»0ll»Pl2»U!2.HHI2»n?»CI?.0 

112 

common/ s  i  npt/  pi.  ui,  Rhi.  of#  xz»  pf.  iri,  i«2 

COMMQn/RAwAV  /  XH( 30) »°P( 30) *CR( 30) »RHR(  30) »EK( 30) »PW(  30) »R 
IMP (JO) 

COMMON  /C1ANU2/  PF  1  $  UF  1  .  RHF  1 .  EF  1 .  C F  1 »  UF  1  ,  PF 2  .  UF 2 .  HHF?  *  FF  ■>  ,  C 
1 F  2  » OF 2  *  XF  ,  X5 
COMMON/llTTS/OIT 
COMMQN/NDIM/TT.  TTMAX.X*Z»TMAX 

cdmmun/ref l/thef  »  r 

COMMOn/GAM/(',AMM(  2  ) 

CUMM(JN/C0N/TC0N#0T2,  PT^L 
CuMMUN/NNEWW/KSrll]CK»KTLLG»  IDT 

905  F  ORM  A  T  (  5E  1 5  •  6  ) 

90/  formahih  ."Initial  the*  o.o"/ 

Uh  *  ”  (  T  )  t  I  mi  of  the  first  time  line  " » 1  PE  1 5 . 8/ 

(lH  »"( TMAX )MAX IMUM  run  TIME  ".1PE15.8) 

906  F(JRMAT(  1H0, "GAMMA  FOR  REGION  UNt "  .  1  PE  1 5 . 8  / 

1  1  H  ."GAMMA  FUR  REGION  I  WO" » 1  PE  1 5 , 6 ) 

9 1 1  K1RMATCIH  . "MEG  1  ON  TwO  PROPERTIES  "/ 

(lH  »6X»"Ul?"» 12X,"CI?"»  1?X,"P12". 11X»"RHI2"» 12X,"EI2"/ 

(  lH  »5( 1PE15.8) ) 

915  F URMA r ( J5A2  ) 

916  F0RMAT( I2.1F15.6) 

91/  f  ormat (  ihi  ."dimensional  input  data"//// 

(lH  ."TIME  units  AHE(".*>A ?»")") 

918  F0RMAT(1H0.S(2X."("«>A?'">">  ) 

919  FORMAT ( 1 H  ."HLG1UN  ONE  PROPERTIES"/ 

(  1H  .6X."UI 1". 1?X."C I  1"* 1?X."PI 1".  1  1X."HHI 1" . 12X."F 1 1"/ 

( 1H  » 5 ( 1  PE  1 5  ♦  6  )  ) 

920  FORMAT  (  1H0."(  I N  )  T  H  E  INITIAL  RAREFACTION  IS  DIVlOEl)  INTO  ".I 

IP.  "  SUn0lV151ONS"/ 

2  1M  »"(  XZ  ^RADIUS  OF  REGION  ONE  " »  1  PE  1 5 , 8 » " ( " . 5 
3A?»")") 

92 1  FuRMAT(lHO."(XZ)HAOIUS  OF  REGION  ONE  ".lPElS.d) 

922  FORMAT!////"  NON-DIMENSIONAL  INPUT  DATA"///) 

925  FORMAT!  1  HO . "  (  PTUL  )  SFFOC*  PRESSURE  JUMP  TOLERANCE  ".1PE1S.8) 
HEAD  90s . T 
OT*T 

HE AO  90S.GAMM( 1 ) ,GAMM( /) 

HEAU  9lS.(  A!  J).J»i.3S) 

011*0.0 

012*0.0 

HEAU  90S. PI  1 .UI 1 . RHI  1 
HEAU  90S.PI2.U12.HHI2 
READ  9l6. IN.XZ.PTUL 
Ell*EuSTEO! 1 .RHI 1 .PI  1  ) 


1 


CllaEuSTCQ(l*EIl*KHIl»r’ll} 

tI2«£Q5TEQ(?»RHl3»PI2) 

CI2aEaSTC0(?»E12»HHI2#Pl2) 

PRINT  9l7,(A( J)»J*1.5) 

PRINT  907.T.TMAX 

PRINT  908»GAMM( 1 ) »GAMM(?) 

PRINT  918»(A< J)»J*6»30i 
print  919»uu»cii.pii.nhii.eii 

PRINT  911,UI2»CI2»PI2#«HI2»EI2 
PRINT  920»IN»XZ.(A< J).U.31.35) 

PRINT  925 » PTUL 

XxZ»XZ 

XZ»1 .0 

TT«T 

DTTaUT 

TTMAX*TMAX 

A ( 1 )«XXZ 

A ( 2 ) ■ XXZ/C 1 1 

A(  3 ) * C  1 1 

A(4)«CU 

A( 5 ) *RHI 1 *C I  1  **2 

A(6)«RHI1 

A<  7)*CI 1**2 

call  nonoim 

PRINT  922 

T*T  T 

uT*ori 

Tmax*ttmax 

PTGL*PTnL*X7 

PRINT  907  »  T  »  TM AX 

PRINT  919.UIl.Cll.PIl.RHl  WEI1 

PRINT  9U  ,UI2»CI2.PI2.PHI2#EI2 

PRINT  9?  1 » XZ 

PRINT  925.PTUL 

PRINT  9?3»(A( J>» J*1 .7) 

923  T0RMAT(//lH  ."CONVERSION  FACTORS  BACK  TO  DIMENSIONAL  UUAnTI 
1TIES"/  1H0. "QUANTITY  MULTIPLY  BY"/ 
t  1M  ."X".9X» 1PE15.8/ 

3  1H  ,"T".8X* 1PE15.8/ 


4 

5 

6 

7  8/ 

OX, 1PE15.8/ 
9,"E".8X»1PE15.8) 
TCUN«A(?) 

CALL  RASH0K1 2. IN. 1 ) 

call  asigpt 

I F  < I0T.EQ.1)UT«0T2/TC0N 
return 


»"U",8X. 1PE15.8/ 

1H  ,"C",8X» 1PE15.8/ 

1H  »"P".8X.  1PL15. 

1H  »  "  R  H  *  •  7 
1h 


end 


subroutine  intfpu  id  i**I3»I4) 

C  SUBROUTINE  To  CALCULATE  PROPERTIES  AT  INTERFACE. 

COHHON/GAIN/  0(2.  1000)  »X(  2#  10U0)  »U(2»  1000)  »C(  2.  1000)  .RH(  2.  1 
1000).  E  ( 2  * 1000) »P( ?  » 1000 ) 

COMMUN/TIMUU/  OT.  UU1#  UU2.  I.  XHU 
100  FQRMATUH  . !  2 . 1  X  .  14  .  3  X  •  6  (  E  1 5  .  »  »  2X  ) . "  *"  .  I  1  .  ?X  »  "  I  NTF"  ) 

110  F URMAT ( 1 H  .10X. "INTERFACE  SOLUTION  DIVERGENT") 

C  1?0  FORMAT!"  *#***"»  6(E1^.S»  2X)»  "*****") 

C2t)<)0  FORMAT!  1H  »  SX  .  "  X  A"  *  1  OX  »  "UA"  .  1  OX  .  "P  A"  .  1  OX  .  "RH  A"  •  VX  .  "E  A"  »  /  , 

C  1 H  . 5E  1 1  * 

C  1  A ) 

C  20n  1  FORMAT! 1H  .SX."XB". 1 0 X • "UB" » 1  OX »" PB" . 1 0 X » "RHH" . 9X . "EH" ./ . 

C  1 H  » 5F  1 1  . 

C  *4) 

xinP!  vi.ov.ox.  Y)*vnnv*y/ox 

L  «  1 

TqL*0.0!)0S 
TOLl  *  l.E-10 
L I M*  30 
K  *  1 

C  *  DEFINE  PROPERTIES  AT  POINTS  1.  2.  3.  AND  4. 


01 

a 

Ut  1  » 

ID 

Cl 

a 

C !  1  • 

ID 

XI 

a 

X!  1  . 

ID 

RH  1 

«  RH!  1 

.  11) 

El 

a 

E!  1  . 

ID 

PI 

a 

P!  1» 

ID 

01 

a 

u!  1  . 

11  ) 

E? 

a 

E!  1  . 

I?) 

P2 

£ 

P!  1  . 

I?) 

RH2 

»  RH!  1 

.  12) 

0? 

8 

U!  1  • 

IP) 

c? 

a 

C !  1 . 

I?) 

U2 

a 

U!  1 . 

I?) 

X  3 

a 

X!  1  * 

13) 

U  3 

a 

U!  1  * 

13) 

C  3 

a 

C!  1 . 

13) 

HH3 

■  RH!  1 

*  13) 

E  3 

E!  1  * 

m 

P  3 

PI  l. 

13) 

X4 

X!  1  • 

I  4  ) 

03 

Ut  1  . 

13) 

U4 

U!  1  . 

14) 

C  4 

Ct  1  . 

I  4  ) 

C  DEFINE  oiff.  quantities 

UX2»X( 1 . I4)-X( 1 . I  3) 
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DU2«U(  1, 14)-U( 1,13) 

DC2  ■  C( 1.  14)  -  C( 1#  *3) 

DE2>E(1»I4)-E(1»I3) 

DRH2*RH(1»I4)"RH(1»I3) 

UP2  ■  P(l,  14)  -  PM,  i  3 ) 

Do2"Q ( 1 , 1 4 )  *0(1,13) 

DXl«Xll,I2)-X(  1,11) 

UUl«U( 1 , I  2  ) -U( 1,11) 

UE1>EM,I2)*E(1,U) 

ORH1  ■  RH(  1 , 12  )  -  RH(lMl) 

DPI  ■  P(l,  12)  -  P<1,  II) 
l)0l«Q(  l,I2)-0(  1,11) 

DC  1  ■  C?  -Cl 

XA  *  (X«A  -  ( U4  -  C4  -  (OUP  -  0C2)*X4/DX2)*DT)/(  1  .  ♦  (HUP  - 
1  DC2 ) *  DT/UX2) 

X 4 A  ■  X3  ♦  U3*DT 

ESTIMATE  pts  a  ano  b 

Xh  ■  ( X4  A  -  cut  ♦  Cl  -  ( OU 1  ♦  DCl)*Xl/DXi)*DT)/(  1.  ♦  (DUW 
1  DC l ) *OT  /OX  1 ) 

DXA  *  XA  -  X3 
DxB  ■  XH  -  Xl 


UA 

X1NPOJ3, 

0U2 , 

DX?» 

OXA) 

CA 

XlNP(C3» 

DC2 , 

0X2* 

DXA) 

PA 

X I N  P  (  P  3  * 

0P2, 

DX?» 

OXA) 

UA 

XlNP(Q3» 

DU2, 

DXP* 

DXA) 

EA 

XlNP(E3» 

0£2, 

0X2* 

OXA) 

hha 

*  EflSTROlL 

,  EA 

,  PA> 

UB 

X 1 NP( u 1 , 

0U1, 

0X1* 

0X8) 

CB 

XInPCCI , 

0C1, 

0X1 » 

OXB) 

PH 

XINP(P1» 

DPI, 

0X1  * 

DXB) 

UB 

XInP(Q1* 

0U1, 

0X1  * 

OXB) 

EB 

XlNPCF.1, 

0E1, 

0X1* 

OXB) 

RHB 

»  EOSTROIL 

,  EB 

,  PHJ 

-  ASSUME  VALUES  FOR  VARIABLES 
U4 A  *  U3 
U46  ■  Up 
C 4 A  «  CA 
RH4A  •  RH3 
E 4 A  *  E3 
H4A  ■  P3 
P4B  ■  P4A 
C4B  «  C8 
RH4B  *  RH2 
E4B  ■  EP 
U34A  ■  Q3 
U24B  ■  02 

1  continue 

BEGINNING  of  ITERATION  LOOP 
X4 A  ■  X3  ♦  U3*0T 
XaB«X4 A 


UMC4A*<U4Afl)A-C4A-CA ) /24 
UPC4b*(U4B*Ud*C4B*CB)/2. 

XA*X4A-UMC4A*0T 

Xd“X4d*UPC4R*DT 

0XA*XA*X3 

L>XB*Xb“Xl 

UA»XInP(U3.()U2.0X2»DXA  ) 

£A*XlNP(E3*n£2.L)X2.nXA) 

PA  ■  XlNP(P3»  0P2.  DX2*  DX A ) 

HHA  ■  E  0  S  T  P  Q 2  #  £  A  »  PA) 

CA  *  EQSTCO( 2  *  LA*  RHA»  PA) 

PEA  «  £<)STPf(2.  E  A*  RHM  ) 

PE  4  A  *  E0STPt(2*  £  4  A  *  "H4A> 

QA*XlNP(0J,nU2.UX2»nXA) 

UH«XINP(U1*DU1»UX1»DXH> 

£H»XlNP(El*n£l*UXl*OXH) 

PH  *  XlNP(Pl#  DPI*  0X1»  DXB ) 

RHB  *  EOSTROC 1  *  EH.  PH) 

Ch  ■  £QSTC0(l»  £H.  RHh*  PH) 

PFtJ  »  Eg  ST  PEC  1 .  Eb.  KHtf) 

P  f  4  b  *  EOSTPt (  1 .  £  4  H .  rtH4B) 

UH*XlNP(Ql»n«l*UXl*OXH) 

04A  ■  03 

04b  »  Q? 

KHC4H  *  ( HHR*CH  ♦  RH4H*C4H)/2. 

RHC4A  a  (HHA*CA  ♦  RH4A*C4A)/2. 

UCX4H  a  {  UR*td/X(J  ♦  U4  U  *  C  4  R  /  X  4  B  )  /  2  , 

UCX4A  ■  (U4*CA/XA  ♦  U4M#€4A/X4A)/2. 

P  RC04B  «  ( PLri*Ob/(  RMH*CB)  ♦  PE4B *G4B / ( RH4R *C 4 b ) ) / ?  . 

P  RCQ4A  a  (P£A.UA/(RHA*CA)  ♦  PE4 A*04A/ ( RH4A*C 4  A ) ) / ? , 

PE04H  ■  (P£B*UH  ♦  PE4H*G4b)/2. 

P£G4 A  ■  ( PE  A  *0  A  ♦  PE4A*Q4A)/2. 

jo  continue 

PRHJ4A  a  (P3/(RH3**2)+P4A/(RH4A**2))/2« 

PRH24d  a  (  P?/(  RH2**2  )  •♦  P4B/ (  RH4b**2)  )/2. 

P4 AP  a  ( Pd/BHC 4b  ♦  PA/RHC4A  ♦  UH  •  UA  ♦  (-(XMU  -  1*)*(0CX4H 

1  ♦  UCX4A)  ♦  P  Rv«04d  ♦  P  RC04A)*DT)/(  1  ./RHC4R  ♦  l./RH 

2  C  4  A  ) 

P 4 B P  *  P4AP 

C  PRINT  2000.XA.UA. PA. RHM.EA 

c  print  ?ooi.xd»ud.PB.RHd.fB 

U4AP  a  (Pd  •  PA  *  RHC  4  0  *UB  ♦  HHC4A*UA  ♦  (-(XMU  -  1,)*UCX4M* 
1RHC4B  ♦  PEB4R  ♦  (XMU  -  1 . ) *UC X4 A*RHC 4 A  -  PE  04 A ) *  0 T  ) / ( H H C 4 
2  A  ♦  RHC4B) 

£4AP  a  E3  ♦  PKH34A  *( KH^A  -  RH  3 )  ♦  Q34A*DT 
E4BP  *  E2  ♦  PRH24d*(  R  H  *♦  B  -  RH2 )  ♦  Q24d*DT 
HH4AP  *  E0STR0(2.  E  4  A  P  '  P4AP) 

HH4BP  *  EOSTROC I .  E4HP  •  P4HP) 

C4 AP  *  E0STCQ(2»  E  4  A  P .  RH4  A  P  »  P4AP) 

C4BP  *  EOSTCOd.  E4HP.  RH4RP.  P4BP) 
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I F  (  AdS((P4AP  -  P4A  )  /P**  A  )  » GT  •  TOL)  GU  TO  3 
1F(  AdS ( U4 A )  .LT.  TOL 1 ^  GO  TO  40 
I F  (  AdS((U4AP  -  U4A  )  /U1*  A  )  .GT.  TOL)  GO  TO  3 

aO  I  F  (  AdS ( (E4aP  -  E4A)/E*A)  .GT.  TOL)  GO  TO  3 

IF(  AtiS((E4HP  -  E48)/E<*8>  .GT.  TUL)  GO  TO  3 

i F (  AdS((«H4AP  •  RH4A)'RH4A)  .GT.  TOL)  GO  TO  3 

I F (  AdS((KH4dP  •  KH4B ) 7  RH4B )  -  TOL)  6’  6»  3 

3  I F ( K  .LT.  LIM)  GO  TO  7 
Print  no 

GO  TO  6 

7  U4  A*  (  U4  AP  +  U4A  ) /2  • 

U  4  0  ■  U  4  A 

C<»A*(C4AP*C4A)/2. 

C4H*(C48P+C48)/2. 

P4A»(P4APtP4A)/2. 

P4B  a  P 4  A 

4  Rh4A*(RH4A.RM4AP)/2. 

NH4H*( RH4H+RH4RP ) /2. 

E4A*( E4A+E4AP)/P. 

L  4  0a  (  E4H  +  E4HP  )  /  2 « 

La? 

PRINT  100.  I»I2»X4A»U4m»C4A*P4A*RH4A*E4A»L 
L  *  1 

PRINT  100.  I.I3.X48.U40.C48.P46.RH4H.E4H.L 

')  K  a  k  ♦  1 
GU  TO  1 

ENU  OF  ITERATION  LOOP 
fa  U(2»l2)aU4AP 
U(2.l3)aU4AP 
P( 2. l?)aP4AP 
P( 2. I 3)aP4AP 
CC  2. I 2)aC4flP 
C( 2# I J)aC4AP 
E( 2. I2)«E40P 
E(2.I3)aE4AP 
H  H ( 2  * 12)*HH4bP 
RH( 2. I  3  )  aRH4 AP 

0(2.  12)  a  0. 

0(2.  13)  a  0, 

X  (  2  .  I  2  )  *  X  A 
X  (  2 . I  3 ) a  X4H 
RETURN 

end 


subroutine  nundim 

common/ ini  T/P  II  .ui  1  .rhA  1  .En.ci1.un.Pi2.u12.RH  12. Ei?.cr>.o 
u? 

CGMMON/NOM/PO.RHO.EO.CO 


CUMMON/nTTSAlTT 

COMhON/NO'iM/1T,TTMAX»X*Z 

HO*Pll 

RH0«RH  1 1 

EO-EI  1 

t 0  ■  C  11 

T  T *  ruco/xxz 

0  T  T ■OTT  *CU/XXZ 

r tmax«ttmax*co/xxz 

P(  1»PI1/(KH0*C0**2) 

KHll  *  RHI1/KHO 

ui i  «  un/co 

til  *  C  1 1  /CO 
1 1  1 *E 1 1 /CO**2 
Pl2*PI2/(RH0*t()**2> 

Hh  I  2  »  R H 1  ? / K H (■ 

U  1  2  ■  U  I  2/CO 
t 12  »  C  I  2/CO 
El?*Ei2/C0**2 
RE  1  URN 
End 


SUBROUTINE  PTARNC 

COMMON /GAIN/  «(2»  1000) »X(  2*1000)  »U(2# 1000) ,C(  2» lOuO) *nh( 2. 1 
1U()0)»  E(2»  1 000  )  #  P (  2* lOUO) 

common  /timuu/  uT #uu l • uu? » I » xmu 

COMMUN/NOCON/  ISl.  I S  2  *  IS3.  ISA,  INTI.  INT2.  IFF 
COMMUN/nTTS/  OTT 
Cijmmon/nnEwr/kshuCk 
DIMENSION  Mf (  1000).  M  A  ^  1000) 

BUO  F0HMATC1H  » I2» IX . 14,  3X'6( E15 ,B»2X ) . 12) 

mil  format  uh  ,iox,i4,"  puints  being  eliminated") 

no 2  FuHHATdH  .10X.I4."  POiNTS  BFINU  ADDED") 

M 1 0  FuKMAKllX*  "  POINT  ELIMINATED  AT  (1#",  U.  ")") 
tt  1  2  FORMAKIIX,  "  NU  PUINT  BEING  ADDED.") 
n 1 4  FUHMATCllX.  "  NEGATIVE  NUMBER  OF  POINTS  ADDED.") 

C  01b  F  uRMA  T ( 1 1 X »  14.  "  P0IN*S  BETWEEN  KK  a".  14,  "  AND  KK  a", 

C  I  A . 

C  1  "  BEING  ADDED.") 

020  F  QRM A  T  C 1 1 X »  "  PUINT  AT  KK  a",  IA,  "  BEING  ELIMINATED  DUE  ID 
I  INTERSECTION  WITH  SHOCK  WAVF") 

02b  FuRMAT  (  IH  . 10X. IA. "POINTS  ADDED  BET.  KK  a  " . 1 A . "  AND  K *  = 
i  "» I  A) 


mo 

Form  A  T  ( ** 

K  15 

c 

040 

FORMATT" 

c 

Hal 

F  nRM A  T ( " 

c 

0  A  2 

F ORKA  T ( " 

c 

Ha  3 

F  OHM A  T ( " 
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C  dsO  FoRMAT("  OTMAX  ■  "#  El*. 8#  "  TE  ■  "#  £15*8) 

0(SO  FuRMAK"  NN  «"#  14#  4  (  E  1  5  »  8  »  2X)) 

1000  FORMAK"  I  SI  ■"»  15#  "  IS?  ■"#15#  "  INTI  -".lb.  "  INT?  =  ", 
1  15#  "  ISJ  ■"#  15#  ■  ISA  ■".  15.  "  IMAX  *"»  I5> 

ClOsO  F  0RMATUE15.8) 

ClOftO  FQRHAK  "  £  ( 1 »  J )  «»#tH,||,  "  PC  1  *  J )  «"#E15.8»"  RH(1#J) 

C  ■"#£15#8# 

C  1  "  X(  1,  J)  «"»E15*8»"  J  ■ "  #  I  *5 ) 

Cltoo  F ORMAT {"  Jl  ■"»  15#  5*.  "  J  ■"#  15.  5X»  "  DT  ■"»  E15.8) 

Cl  1 10  FuRMAT{"  OT  «"#  £15*81 

C1200  FoRMATC"  IFF  «"#  15#  1  OX »  "  OT  ■"#  E15.8) 

C1JO0  FORMATt"  C  (  1  #  J )  ■"#  £l5*8#  "  E(i#J)  ■"#E15.8»"  P(1#JJ  =”• 

C  E  1 5  *  8  # 

C  1  "  RH( 1 #  J )  ■"#£15*8#"  X( 1 » J)  *M#E15*8#n  J  *"#I5) 

XLINP(V1#  OV#  OX#  DY)  *  VI  ♦  DV*DY/OX 

0 1 N  P  (  V?#  nvl.  U  V2#  0  Y 1  ■  M2  *  (DV2*DX1**2*D\/1*DX?**?)*  liY 
<  0X1  *DX2*(  OX  1  *  0X2))  ♦  CDVl'OX?  ♦  DV2*DX1)*(  iiy 

*  )**2/(nxi*nx2  *(oxi  ♦  ox2)) 

c  I  F  (  1  .  G  E  •  3>  0  T  1  *  U  T  T 

1FU  #£q,1)  0T1*DT 
1FC0T  « i»T  «  0  T  1 )  DTl«f)T 
C  OTT  «  .005 

C  OT 1  ■  0T*6*/5. 

C  I F ( I  .Eu.  1)  OT 1  *  OT 

C  1 F ( 0 T 1  .GE*  OTT)  DT  1  «  OTT 

T  MIN  ■  0.8 

tmax  ■  ?.o 
TmAX  ■  2.0 
thin  ■  0,5 
Tmax  ■  100.0 

T MAX*5«  0 
Thin  «  0.9 
tmax  ■  3.0 

0TM1N  ■  TMIN*0T1 
UTMAX  «  TMAX*0T1 

ioi  format ( 1 h  »"dtmin*".ei*,8#"dtmax»".ei5*8) 

1FF1  .  IFF 
Ik  ■  0 

IS1M1  «  isi  -  1 
KKM1  ■  1 
1S1M2  «  ISI  -  2 
lF(KSHOCK.Eo.l)  GO  TO  <<00 
IF(IS1M2  .LE.  11  IS1P  *  ISI 
IFUS1M?  .Lf.  1)  GO  TO  15 
DO  10  KK  ■  ?#  IS1M1 
1F(  IK  .EU.  (  ISI  •  3) )  GO  TO  1 7 
9  TE  ■  ADTCKK#  KKM1) 

1 F C  T E  *GT.  OTMIN)  GO  TO  8 

Ik  ■  ik  +  l 

M£( IK  1  a  KK 


GO  TO  10 

8  KKM1  a  KK 

10  C  ON  T  I  NUf 

IS1P  a  IS1  -  IK 

IF  (IK  ,E0.  0)  UG  TO  U 

IF(MEUK)  .EU.  1 S 1 H 1  >  un  TO  15 

1  1  Tf  *  ADT( 1 S 1  »  1 6 1 M 1  ) 

IF(TE  « fiT •  D TM I N )  GO  TO  15 
IK  «  IK  ♦  1 
ME  C  I  K  )  *  I  SI  Ml 

17  CONTINUE 

1S1P  a  IS1  -  IK 

15  CONTINUE 
IS2P1  *  I S2  ♦  1 
INT1M1  .  INTI  -  1 

I F (  l  S2P 1  ,E0.  I N 1 1 )  GO  TO  24 

KKM;  a  I S2 

G(J  Tl?  401 
400  I S2P 1 *2 

InTIMUinTI-1 

ISlP«Ibl 

401  continue 

UQ  20  KK  ■  I S2P 1 »  I  NT  1 M 1 

16  TE  «  AO  T( KK  #  KKMl) 

IF(TE  «r,T ,  n Tm i n )  go  to  iy 
Ik  *  IK  ♦  i 
ME(IK)  a  K K 
GO  TO  20 

iy  KKMl  a  KK 
2U  CONTINUE 

22  continue 

InTIH  *  InTI  -  IK 

IF  l  IK  « E  Q «  0}  GO  TO  *3 

1FCMEC IK)  ,f«.  INT1M1)  GO  TO  25 

2  3  Te  *  ADT ( Inti*  intimi  ) 

I F (  T E  *GT.  nTMlN)  go  TU  25 
IK  ■  1K  ♦  1 
me ( l k )  a  intimi 

2 4  InTIP  a  INTI  -  IK 

25  CONTINUE 

I Nl 2P 1  a  I  NT  2  ♦  1 

IS3M1  a  IS3  -  1 

IFCINT2P1  * EO •  I$3>  GO  TO  32 

KKMl  a  InT? 

IF  (  INT2P1.EQ.  IS3MI  >  GO  TU  31 
UO  30  KK  »  INT2P1#  I S  3 M 1 
?6  TE  *  AOT ( KK »  KKMl) 

I F (  TE  .GT.  OTMIN)  GO  TO  ?9 
Ik  *  IK  ♦  1 
ME ( l K  )  a  KK 


GO  TG  30 

V9  kkmi  ■  kk 

30  continue 

31  continue 

IS3P  -  I S 3  -  IK 

IF  <  IK  .EQ.  0  )  QO  TO  34 

IFCME(Ik)  .Efl.  1 S 3 H i )  «n  TO  35 

34  TE  ■  ADT(  I  S3.  IS3M1) 

IFCTE  • GT  .  0  T  M I N }  GO  TO  35 
IK  *  IK  ♦  1 

M£(lK)  >  IS3H1 

32  IS3P  ■  I S 3  -  IK 

35  continue 

IS4P1  *  IS4  ♦  1 
I FFm 1  *  IFF  -  1 
KKH1  -  I S  4 

00  40  KK  *  IS4P1.  IFFMl 

36  TE  ■  ADT ( KK  »  KKMI) 

IF(TE  *  GT  •  DTmIn)  GO  TO  39 
IK  «  IK  ♦  1 
ME(lK)  «  KK 
GO  TO  40 

39  KKH1  *  KK 

«o  continue 

IFFP  .  IFF  -  IK 
IF( IK  .EO.  0)  GO  TO  41 
IF( ME( IK )  *EU.  IFFMl)  uq  tO  45 
4i  TE  a  AOT( IFF  *  IFFMl) 

IF(  r£  «GT,  oTMIN)  GO  TO  45 
IK  *  IK  ♦  1 
ME(IK)  a  IFFMl 
IFFP  a  IFF  -  IK 

45  continue 

1  S  1  ■  IS1P 

c  -  END  OF  ELIMINATING  POINTS. 

I S2  a  I S 1 P  ♦  1 
I S  3  a  IS3P 
I S 4  a  I S  3P  ♦  1 

1nT2  *  intip  ♦  1 
1 N  T 1  «  INTIP 
IFF  a  IFFP 

IF  C  IK  .EQ.  0)  GO  TO  56 
JK  a  1 
K  a  1 

00  55  KK  a  1.IFF1 

IF  (  ME(JK)  .EO.  KK  )  'JO  TO  53 

P(  1.  K  )  aP(  1  .  KK) 

U  (  1 .K )  a  U( 1 .KK) 

KH(  1  .  K  )  ■  KH(  l.KK  ) 
t( 1 .K)  a  E(  1  .KK) 
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Cl  1 »  K  )  ■  Cl  1 »  KK  ) 

Ul  l.K )  a  Qt 1  »  KK  ) 

XI  l.K)  ■  XI 1  »  KK  ) 

K  •  K  ♦  1 
GO  TO  5S 

si  Continue 

IF  l  JK  ,E().  IK)  GO  TO  S') 

JK  *  JK  ♦  1 

55  CONTINUE 

56  CONTINUE 

1FIIK  .EO.  0)  GU  TO  49 
00  4tt  J  *  1 .  IK 
4tt  CONTINUE 
uv  CONTINUE 
C  1)1)  SO  J  ■  1»  IFF 

C  SU  PRINT  dOO  »  I  » J *xl 1 » J)  .ul  1 » J) »Cl 1  * J)  »Pl  1  » J)  »RHl 1 » J) »El  1  »  J  )  * L 
C  -  AOO  POINTS  STARTS. 

C  -  IK  ■  TOTAL  NUMttEH  OF  POINTS  TU  dE  A U D E 0 * 

C  I F t IFF  ,GT.  0  )  GO  TO  100 

l  K  *  0 

IFIKSHOCK.EQ.  1  )  GU  TO  *?0 
OQ  6 U  KK  *  ?»  151 

KkMI  *  KK  -  1 
MA(KKMl)  a  IK 
TE  *  AOTIKK.  KKMI) 

c  print  8so»i)Tmax.  te 

I  F l X l l #  KK )  .LT.  0.9E0**t 1 .IFF)  )  GU  TO  60 

lFCXII.KK)  .LT.0.99E0*A(1,IFF))  GU  TO  60 

UTMAX*3.OE0*OT 

IFI TE  .LT .  OTMAX)  GO  TU  60 

NN  >  l  X  l  1  #  KK)  -  X  (  1  •  AKMl))/UCtl»  KK)  ♦  C  l  1  *  KKM1))*0T)*? 
1  . 

IFINN  »LE.  1 )  GU  TO  60 
IK*  IK  ♦  NN  •  1 

ko  continue 

MAC 1S1  )  ■  IK 
1S1P  *  I  SI  ♦  IK 
IS2P1  ■  IS?  ♦  1 
GO  TU  4?1 
4  ?()  I S2P 1 *2 

4?l  continue 

00  65  KK  *  IS2P1.  INTI 
K  K  M 1  *  KK  -  1 
MAlKKMl)  *  IK 
TE  ■  AOTIKK.  KKMI) 

OTMAX* 1 . 8E0*UT 

IFITE  .LT.  OTMAX)  GO  TU  65 

NN  *  t  X(  1  .  KK)  -  XU.  KKMl))/l  ICl  1.  KK)  ♦  C  l  1  .  KKM  1  )  )  *1)  T  )  *  ? 
1  . 

IFINN  .LE .  1  )  GU  TO  65 
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C(  1  »  KKMl  ))*0T  )  *'? 


IK  ■  IK  ♦  NN  •  1 

(s*>  continue 

na( inti )  a  ik 

intip  ■  InT 1  ♦  IK 

1NT2P1  •  INT2  ♦  1 

00  70  KK  ■  I N  T  2P  1 »  I  S 3 

KKMl  a  KK  -  1 

MA(KKMl)  ■  IK 

TE  ■  AOT(KK»  KKMl  ) 

!)TMAX*3«0E0*UT 

1F( TE  # LT .  DTMAX)  GO  TO  70 

NN  ■  (X(l,  KK)  -  X(l,  KkMI ) )/( ( C( 1 »  KK)  ♦  C(l»  KKMl)) 
1. 

)■?  PQRMATdH  ,"KK*"»  l4,"NNa",l4»"Tta",E15.8»,’DTa",E15*8) 
I F  C  NN  .LE.  1)  GU  TO  70 
IK  «  IK  ♦  NN  •  1 

70  continue 

M  A ( I S  3 )  a  IK 

1S3P  ■  I S  3  ♦  IK 

IS4P1  a  IS4  ♦  1 

00  75  KK  *  I S 4 P  1  *  IFF 

KKMl  «  KK  ■  1 

MA(KKMl)  a  IK 

I F  C  IFF  ,gT.  0.)  GO  TO  ^5 

TE  a  ADTCKK,  KKMl) 

IFCTE  .LT.  OTMAX)  GO  TU  75 

NN  a  <  X «  1  .  KK)  -  X<1.  KkM1))/((C(1»  KK)  ♦  C(l»  KKMl)) 
1  • 

1 F ( NN  .LE.  1)  GU  TO  75 
IK  *  IK  ♦  NN  •  1 

75  continue 

ma(IFF)  a  Ik 
1FFP  a  IFF  ♦  IK 

) QRM A  T (  1  H  #**lFFa»,  l4#”lFFPa",  I4»"!Ka**»  14) 

1F(MA( IFF)  .GE .  1 )  GO  ]0  76 


TU  75 

KkM  1 ) ) / ( ( C  (  1  » 


C(  1,  KKM 1  )  ) *D  r  )*? 


GO  T 0  77 
X(l.  1FFP) 


X(  1  » 

P  (  1 »  IFFP)  a  PCI*  IFF) 

U( 1 #  IFFP)  a  U( 1 .  IFF  ) 

L { 1 »  IFFP)  a  E ( 1 •  IFF) 
RH(1#  IFFP)  *  RH  (  1  »  I  F »“  ) 

C  c  1  a  IFFP)  a  C( 1 f  IFF) 

UC  1  *  IFFP)  a  0(1,  IFF  ) 
CONTINUE 
K  a  IFFP 

on  60  KK  *  ?  *  I F F 

IN  a  IFF  “  KK  ♦  1 

NFT  a  MA(  IN  ♦  l )  -  MAC  IN) 

K  a  K  -  NPT  -  1 

X(  1,  K)  a  X(  l ,  IN) 


IFF) 

IFF) 

IFF) 

IFF) 


M(  1,  K)  a  P(l,  IN) 

<J(  1,  K)  a  UC  I,  IN) 

E(  1»  K)  *  £(  1  ,  IN) 

KH(1,  K)  a  HM( 1 ,  IN) 

C(i.  K)  a  cc  1  ,  IN) 

U(  1  »  K )  a  0(1,  IN) 

ho  continue 

CALCULATE  PRUPEHTIES  At  ADDED  PlS. 
KCHCK  a  i 

IE< KSHOCK.Ey.  1 )  KCHC**2 
K  ■  1 

7V  0(1  TO  (Hi,  H?»  03,  HA  )  »  KCHCK 
Hi  HI  «  2 

*2  *  ISl  -1 
GO  TU  VO 
H?  Ml  a  IS?  ♦  1 

If(  K6H0CK ,FU.  1  )  Ml*? 

M2  *  INTI  -1 

IECM1  ,EQ.  INTI)  GO  TO  92 
GO  UJ  VC 

MJ  Ml  a  1 N  T  2  ♦  1 
M?  a  IS3  “1 

I F  C  H  3  ,E0.  I  S3)  GU  TO  V} 

GO  TO  VO 
H 9  Ml  -  ISA  ♦  1 
M2  ■  IFF  -1 

IF(M1  .fo.  IFF)  GO  TD  Vfl 
90  CONTINUE 

UO  06  KK  a  Ml,  M2 

NPT2  a  M  A ( K  K  ♦  1)  -  MA 1 KK ) 

NPT1  a  MAC  KK  )  -  MAC  KK  "  1 ) 

K  a  K  ♦  NPTl  ♦  1 
I  F ( NP  T  1  ,EQ.  0)  GO  TO  Oft 
IFCNPT1  ,LT.  0)  PHINI  0 1 A 
XP1  a  NPTl  ♦  1 

UX  a  X  C 1  *  K)  -  X( 1,  K  *  NPTl  -  1 ) 

0  Y  1  »  “DX/XPl 

Uxl  a  Ox 

UX2  a  X (  1 ,  K  ♦  NP  T  2  ♦  * )  -  XC  1  ,  K) 
OY  a  0  Y  1 

00  05  IA  a  l,  NPI 1 
^2  »  PCI,  K) 

U2  ■  U( 1 ,  K) 

E2  «  L ( 1 ,  K ) 

MH2  a  Rh( 1 ,  K ) 

02  a  g i i ,  k  ) 

C 2  •  CC  l ,  K) 

UP1  =  P( 1 ,  K)  “  PCI,  K  •  NPTl  “  1 ) 

DU  1  -  U C 1 ,  K)  -  U(l,  K  -  NPTl  •  1) 

DEI  *  F.Cl,  K)  -  E(l.  K  -  NPTl  -  1) 
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L) R H 1  «  RH  (  1  »  K)  -  RH(1»  K  -  NPT1  -  1) 

OQl  ■  Q( 1.  K)  -  QM»  K  -  NPTl  -  1) 

UC 1  ■  C ( 1  *  K)  -  C(l#  K  -  NPT 1  -  1 ) 

DP2  ■  P  (1  #  K  ♦  NPT2  ♦  1  )  -  PM#  K) 

L)U2  ■  U(l»  K  ♦  NPT2  ♦  1)  -  U(l»  K) 

0£2  ■  E(l>  K  ♦  NPT 2  ♦  1)  -  EM#  K) 

URH2  ■  RH (  1  *  K  ♦  NPT2  *  l)  -  RH(1*  K) 

0Q2  ■  0(1*  K  ♦  NPT2  ♦  1)  -  Q(  1  *  K) 

UC2  *  C(l»  K  ♦  NPT2  ♦  M  -  C(l*  K) 

PCI#  K  -  IA)  ■  UINPCP2*  0P1*  0P2.  DY) 

U(l*  K  -  IA)  *  U I NP ( U2  *  nut*  UU2#  OY) 

t  (  1  *  K  -  IA)  ■  U I NP ( £2  *  0£t»  DE2#  OY) 

K H ( 1 •  K  -  IA)  ■  QINPCR"?#  ORHl*  DRH2*  DY) 

U(l#  K  -  IA)  ■  UlNP(02»  001*  002 »  DY) 

1(1#  K  -  IA)  ■  UINPCC2*  DC  1 •  UC2#  DY) 

X( 1#  K  -  IA)  ■  X( 1 »  K)  ♦  DY 
X I A  *  IA  ♦  1 
l)Y  «  DY  1  *X I  A 

PRINT  800#  I»J#X(  1# J).U(1,J),C(  i»J)»P(  l#J)»RH( l.J).E(  1# J)#L 
*5  CUNTlNUE 
J  «  K  -  IA 
H6  CONTINUE 

GO  TD  (01#  02#  93#  9A)»  KCHCK 
Ol  M3  *  IS1 
GO  TO  96 
g2  M3  «  INTI 
GO  TO  96 
oJ  M 3  *  I S 3 
GO  TU  96 
oA  GO  To  9fl 
06  CONTINUE 

NPT  2  ■  MACM3)  -  MACM3  “  1) 

K  ■  K  ♦  NPT2  *  1 
1FCNPT2  .EE.  0)  GU  TO  *7 
XP2  *  NP T 2  ♦  1 

Ox  ■  X( 1 #  K)  -  X  C  1  #  K  “  NPT  2  -  1  ) 

UYl  *“Ox/XP? 

OY  =  UYl 

00  95  IA  «  1#  NPT2 
P 2  *  PCI#  K ) 

U2  ■  U( 1 .  K  ) 

E?  *  EC  1 »  K) 

H H 2  *  RHC 1 »  K  ) 

02  *  (J (  1  •  K) 


C2  =  C(l# 

K) 

0P2  *  PC  1  # 

K) 

m 

PC  1  # 

K  -  NPT2  “ 

1  ) 

002  ■  U(  1  . 

K  ) 

m 

UC  1  # 

K  “  NPT?  • 

1  ) 

0£2  ■  E(  1# 

K) 

m 

E(1  » 

K  •  NPT2  - 

1  ) 

0RH2  *  RHCt. 

K  ) 

-  RHC  1 »  K  -  NPT2 

m 

UQ?  «  «(  1  . 

K) 

- 

UC  1  > 

K  -  NPT?  * 

1 ) 

*  ! 


UC2  ■  C ( 1 •  K)  -  C(l»  K  -  NPT2  -  1) 

J  ■  K  -  IA 

P(l,  J  >  ■  XLINP(P2»  DP2.  DX*  DY) 

U(  1  *  J  )  ■  XL  I NP ( U2  *  0U2 •  OX#  DY) 

E(l#  j  )  -  XLINPCE2*  0E2*  DX*  DY) 

U(l#  j  )  ■  XLInP(Q2»  DU2*  DX*  DY) 

C(l»  J  )  *  XLINPC C2»  DC2*  DX*  DY) 

HH( 1  *  J  )  ■  XLINP(NH2»  0RH2*  DX*  DY) 

X(  1*  J  )  ■  X( i»  K)  ♦  DY 

X I A  *  IA  ♦  1 
UY  s  D Y 1  *  X  I  A 

C  PplNT  800*  I  # J* X( 1 # J) *u( 1 » J) »C(  1  » J)  *P(  1  * J) *RH(  1 » J) >E(  1  * J) *L 

s»5  continue 
g 7  cqntinuc 

k  ■  K  ♦  1 

kchck  ■  KCHCK  ♦  1 

1HKCHCK  .LE.  4 )  GO  TO  79 

c  -  complete  ado  pts. 

gti  CONTINUE 

00  101  J  a  ?»  IFF 
! F( MA( J )  .EO.  MA(J-l))  GO  TO  101 
N  ■  M A ( j )  -  MACJ-1) 

JMl  *  J  *  1 

101  continue 

i  s  i  «  i  s  1  p 

1 52  *  ISIP  ♦  1 

1 5 3  ■  IS3P 

154  *  I S 3P  ♦  1 
lNl'l  «  I  NT  1  P 
INT2  «  INTIP  ♦  1 
IFF  *  IFFP 

C  PRINT  lOOO#  ISl*  IS2.  1  NT  1 »  INT2*  IS3*  IS4*  IFF 

(  on  99  J  >  1.  IFF 

C  g9  PRINT  800*  I  * J » X ( 1 « J ) » U( l , J ) » C ( 1  * J ) » P l 1  * J )  * HH( 1  * J ) * E ( 1  * J ) • L 

loo  continue 

C  on  250  J  a  1*  IFF 

(.  PsO  PHINT  800#  I  »  J  #  X  (  1#J)*U(1#J)»C(  1  *  J  )  *  P  (  1»J)*RH(1*J)»E(  1#J)»L 
X  I S 1  ■  X(  1 .  ISl  )  ♦  UU 1 *DT 1 
XISJ  ■  X(l.  1 S  3 )  ♦  UU?*0T1 
X ( S2  a  XIS1 
XIS4  ■  X  I S  3 

XINT  ■  X(l»  InTI)  ♦  U(l#  INTI )*0T1 
C  “  IK  ■  NO  PTS  ELIMINATLO. 

C  -  ME ( IK )  ■  PT  OF  URIGI^AL  NO  BEING  ELIMINATED 

Ik  ■  o 

I F C KSHOCK.EO*  1  )  GO  TO  116 
IS1M1  a  ISl  -  1 
IF  (  UU1  .GT.  0.  )  GO  TO  107 
OQ  105  KK  ■  2*  1S1M1 

XD  ■  XlSl  -  (X(l»  KK)  +  U(l»  KK)«UT1) 
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1 F ( X  D  « OT  «  0*)  GO  TO  IO5 
IK  »  IK  ♦  1 
M£ ( 1 K )  a  KK 
105  CONTINUE 

IS1P  ■  IS1  -  IK 
GO  TO  116 

107  IS1H  ■  IS1  •  IK 
110  CONTINUE 

1S2P1  •  IS2  ♦  1 

1NT1M1  «  INTI  -  1 

IF  (  UU1  .LT.  0.  )  GO  <0  116 

UO  115  KK  •  1S2P1.  I  NT  1 M 1 

XO  ■  XlSl  -  (  X  ( 1  »  K  K  )  ♦  U ( 1  *  KK )  *  DTI  ) 

X[)  *  XlSl  -  ( X ( 1 »  KK  )  ♦  U  (  1  »  KK  )  *  DTI  )  -0.001D0 
IF(XD  .LT.  0.)  GO  TO  US 
Ik  *  ik  ♦  l 
M£(IK)  *  KK 

1 1 5  continue 

116  CONTINUE 

1NT1P  *  INTI  -  IK 
1?0  CONTINUE 

INT2P1  «  I NT2  ♦  1 

IS3M1  *  I S  3  -  1 

IF  l  UU2  .GT.  0.  )  GO  ^0  127 

00  125  KK  *  1NT2P1#  ISJMt 

XD  »  X  IS  3  -  l  X  (  1  • K  K )  ♦  U( 1  *  KK )  *  DTI) 

XI)  ■  X  I  S  3  -  (  X  (  1  »  K  K  )  *  U  (  1  #  KK  )  *  DTI)  -0.00100 
I F ( XD  .GT.  0.)  GO  TO  12s 
Ik  *  IK  ♦  l 
M£CIK)  a  KK 

i ?5  continue 

I  S  3  P  a  I S  3  -  IK 
GO  TO  136 

\o7  IS3P  a  ISi  -IK 

no  continue 

IS4P1  a  1 S  4  ♦  1 

IFFM1  a  IFF  ••  1 

IF  (  UU?  .  LT  »  0.  )  GO  >0  136 

00  135  KK  »  IS4P1  #  I F  F  M 1 

XD  «  X I S4  -  (X(1#KK)  ♦  U(h  KK)  *  DTI) 

l  F  ( XU  .LT.  0.)  GO  TO  13«> 

lK  *  IK  ♦  1 

Me  <  IK)  a  KK 

i  is  continue 
3ft  continue 

1FFP  a  IFF  -  IK 

i 40  continue 

IF  (  IK  « EG .  0)  GO  TO  Iftl 


o  o  o  o 


DO  150  KK  *  1.  IFF 

IF  l  ME(JK)  .EG.  KK)  GU  TQ  145 

P<  1  »  K  )  ■  P( 1 »  KK  ) 

U(  1.  K)  •  U(  1.  KK) 

HH(1»  K)  ■  RH(  l»  KK) 

E  (  1 »  K)  ■  E  ( 1 »  KK) 

C(  1  »  K)  »  C ( 1 »  KK  ) 

0(1*  K)  ■  0(1#  KK) 

X< 1 .  K)  ■  X(  1  »  KK  ) 

K  «  K  ♦  1 
GU  10  150 

U5  continue 

I F  ( JK  «EG<  IK)  GO  TO  00 
JK  ■  JK  ♦  1 
150  CONTINUE 

lFUK  »F.  Q  •  0)  GU  TO  1 6  i 
1 a6  CONTINUE 

00  160  J  *  1»  IK 
1*0  CONTINUE 
1 6 1  151  *  IS1P 

152  ■  IS1P  ♦  1 

1 5 3  «  1S3P 

154  ■  IS3P  +  1 
INTI  *  I  NT  1 P 
1nT2  «  I  NT  1 P  ♦  1 
IFF  *  IFFP 

C  UQ  251  J  »  1  »  IFF 

(  25 1  PRINT  800*  I  > J » X (  l » J ) » U< 1 . J ) » C ( 1 » J ) »P(  1 • J ) » R H ( 1 » J ) # E ( 1 • J ) * L 

return 

lno 


SUbROUT  I  NE  RASHOKlL.  I  *>»  *  HSR) 

COMMON/GAIN/  U (  2  *  1 000 ) »  X  ( 2 » 1 000 ) » U( 2 » 1 000 ) »  C  ( 2»  1000) » KH  ( ? » 1 
1000).  E<?» 1000) »P< ?» 1000) 

C0MMUN/INIT/PI1  .UIl#HHll.EIl.CIl.QIl»PI2*UI?.RHI2.EI2.C0.O 
II? 

common  /  shk4o  /  u4r»c**h»hh4B»E4B»P4B»X4B»q4b  »uuu 
CuMMUN  /SHK4A  /  U4A»C4*.RH4A»E4A»P4A»X4A»U4A 
CllMMON/TlMUU/  DT,  UUi.  UU2.  1#  XMU 

CQMMON/RAwAV  /  XR( 30)  »UR( 30) »CRC 30) »RHR(  30) »FR( 30 ) »PR( 30 ) .R 
1  HP  I  30  ) 

COMMON  /C1AN02/  PF 1 . UF I . RHF 1 . EF 1 . CF 1 . GF 1 # PF2 . UF? » RMF2 . I F ? . C 
1 F  ? » Q  F  2 »  XF.XS 

COMMON/SINPT/  PI.  UI.  HmI.  UF»  XZ.  PF»  IR1.  IR2 
IN  -  *  OF  INTERVAL*  OIV.  for  PRESSURE..  MSR  ■  1 

RIGHT  S 

and  left  rawave.  msr  ■  2  -  left  shock  anu  right 

RAWAVE. 
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N I  T ■  1  00 

TOL  ■  0,000s 
hax«ioo 
1 R  l  «  1 
I R 2  *  IN  t  1 
K  ■  1 
PI  •  PI1 
UI  ■  U I  1 
P4A  *  P I  2 
HHi  *  RHII 
NH«A  «  RHI? 

L4A  ■  El? 

C4A  ■  C 1 2 
U4 A  *  UI? 

MS  *  1 
MK  *  1 
LA*-1 
LS-2 

IF  (MSR  « El) ,  1)  GO  TO  i0 
MS  *  2 
MR  ■  2 
LA*? 

LS*1 

10  CONTINUE 

PF  *  (PI  ♦  PA  A ) / 2 , 

P 4 B  •  PF 

CALL  HAwAy/E  (LA, IN, MR) 

CALL  SHOKLQ( LS.MS,  1  ) 

\b  CONTINUE 

U4H  ■  (UAB  ♦  UF)/2, 

call  smoked  (ls»ms,2) 

PF  ■  PAH 

CALL  RAwAVE  (LA, IN, MR) 

IF  ( AbS(  (UAH“UF)/UF)  ,*-T.  TOL)  GO  TO  30 
IF  (K  ,LT.  NIT)  GO  TO  ?0 
GO  TO  30 
CONTINUE 
GO  TO  IS 
iO  CONTINUE 
INF  ■  If)  ♦  1 
00  AO  J  *  1,  INF 
X( 1 , J)  >  XR<  J) 

P( 1, J)  s  PR( J) 

U( 1, J)  *  UR( J  ) 

L ( 1, J)  *  ER( J) 

RH(  I  ,  J)  ■  RmR( J ) 

C( 1 , J)  *  CR( J) 

U( 1 , J)  r  0, 

4O  CONTINUE 

l NT  1  *  INF  ♦  1 
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INT2  ■  INTI  ♦  1 

X( 1 « InT  1  )  a  XZ  ♦  UF  *0T 

P(  1  »lN7l  )  •  PF 

NH( 1  *  INTI  )  •  HHHC  INF  ) 

E (  1  *  I  NT  1  )  ■  t H ( INF) 

U(  1»InT1  )  «  UF 
CC  1#INT1 )  a  CK( INF) 

U(  1#  I  NT  1  )  •  0. 

X( 1 * InT2  )  •  X( 1 # INT  1 ) 

P( 1 » InT?)  ■  PF 
HH(1#INT2)  *  HH4b 
L (  1* I N T 2  )  a  E4B 
U( 1* INT  2)  a  UF 
C( 1# INT?)  a  C48 
U(  1  * INT?)  a  0. 

1  S  3  a  I  N  T  2  ♦  1 
IS4  *  I S 3  ♦  t 
X(1*I53)  a  xZ  ♦  UUU*0T 
P(  1* 1 b  3  )  a  PF 
KH( 1  *  I  S3 )  »  PH4d 
U(  1,1  S3)  «  UF 
E( 1* IS3)  ■  Md 
U( 1 » IS3)  «  0 . 

C( 1. 153)  a  C4d 
X( 1# I S4  )  a  Xl 1 .  I  S  3  ) 

P(  1*  ISA)  a  H>4A 
KH(1»IS«)  a  HH4A 
E ( 1  *  I  54  )  a  £<*A 
C ( 1 » 154)  a  C4A 
U( 1 * I  54 )  a  U4A 
U( 1* 154)  a  0. 

I  X  *  1 

X  F  *  X  <  1  *  InTI) 

<S  ■  X(  1  *  I  S 3  ) 

PF i  *  P(  1  »  INTI) 

UF)  «  U l  1  *  INTI) 

NhFI  a  R  H  C l •  INTI) 

EF1  ■  E(  1.  INTI) 

CF1  *  C(  1  .  INTI  ) 

QFi  *  0(1*  INTI) 

PF2  a  P(l,  1 NT2 ) 

2  a  U(  1 .  I  N  T  2 ) 

HHF2  a  RH(  l.  I  NT  2 ) 

EF2  a  Ell*  I N  T  2 ) 

CF2  «  C( 1  •  INT2) 

OF  2  ■  0(1*  I  NT  2 ) 

V 1 4  FnRHAT(lH0.I5»SX.6( lPElSifl)) 

9 ?2  FORMAT!  lHl. "THESE  POlN^  DEFINE  THE  INITIAL  SINGULARITY  AMD 

1UT  THE  INTERFACE”/ 

2  1HO»"POINT  NO.  "t 7X»"X" 0 14X»"U" » l4X*"C" *  1 4 Y," 
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JP"»  1  3X."RH"»  UX."E"/  lf1  ."RAREFACT  ION  IN  REGION  ONE") 
FnRMAKlH  ."INTERFACE") 

92b  F'URMATdH  *  "SHOCK  IN  REGION  TWO") 

*26  F ORHAT( 1HO. "SHOCK  VELOCITY  ".1PE15.8) 

PRINT  V?? 

A*  1  «  0 

DO  50  J« 1  *  I S A 

IF( J.EQ.lNTl JPRINT  92A 

IF(J#E0»IS3) PRINT  925 

PRINT  914»J»A»U(1»J).CI1»J)#P(1»J)»RH(1.J)#E(1»J) 

SO  CONTINUE 

Print  926. uuu 
return 

E  NO 


SUBROUTINE  RAwAVEIL.  IN,  HR) 

c  -  compute  rarefaction  WAVE. 

(  -  RHI,  PI.  on  INITIAL  properties. 

c  ■  p m  pressure  at  contact  line. 

c  -  INI  NO.  OF  INTERVALS  TO  HE  DIVIDED  FOR  PRESSURF. 

(  -  MR  *  1,  LEFT  RUNNING  WAVE)  MR  *  2.  RIGHT  RUNNING  WAVE. 

CuPMQN/RAwAV  /  XI  JO  )  .U<  30)  »C(  30) .  RHI  30)  »E(  30  )  »P(  30  ) »  RHPMO  ) 
COMMON/SINPT/  F»»  01.  "Hi.  UF *  XZ»  PF »  I R 1 •  IR2 
COMMUN/T I MUU/  OT.  UU1.  UU2 .  I#  XMU 
100  FORMATl"  NOT  CONVERGING  IN  RARAVE  ON  N  ■" »  12) 

V  oO  FfjRHAT(l2x.  "X".  1 6  X .  MU".  16X.  "C".  I7X.  "P".  17X.  "HH".lf 
IX,  "L"  ) 

lunO  EuRMATUx.  ft  I  E  1  b  .  rt  »  ?X).  I?) 

1  (>  1  F  0  P  M  A  f  (  1  H  ."ROUTINE  RAWftVE  L*"»I3»"  IN*". 13."  MR*". 13) 

F  OL  *  0,0001 

MAX* 1 00 
SIGN  *  - 1  , 

IF(MH  ,F0.  ?)  SIGN  a  1, 

r;  -  l)t  F  INf  PROPERT  IES, 

P(  1  )  •  PI 
U( 1 )  ■  ul 
RHI 1 )  =  RHI 

ECl)  *  EOSTFUI  L.  RHI  .  ,J  I  ) 

C( 1 )  *  EOSTCUIL.  E (  1  ) .  RHI .  PI ) 

X ( 1 )  a  XZ  ♦  IUI1)  ♦  SIUN*C(1))*0T 
x i n  *  In 

OP  *  (PF  -  Pi ) / X  1  N 
I N  F  »  IN  ♦  1 
00  5  N  a  2,  INF 
b  P ( N )  *  P(N-l)  ♦  OP 
XINF  a  InF 

00  25  N  ■  ?,  INF 

C  -  ASSUMf  VALUES  FUR  VARIABLES. 


on 


K  ■  1 

RH! N)aRH! N-l )  *  (  1  .♦ 1*/XINF*SIGN) 

C(N)  a  C(N-1 ) 

KhC  ■  RH( N-l )*C(N-1 ) 

c  -  calculation  ok  properties* 
io  continue 

NHC  •  (RH( N-l )*C(N-1 )  ♦  RH(N)*C(N))/2> 

PRH  a  (P(N-1)/CHH(N-1)**2)  ♦  P(N)/(RH(N)**?))/2* 

U  ( N )  a  SlGN*(P(N)  -  PP-lll/RHC  ♦  U(N-l) 

E  (N)  a  PRH*!RH!N)  -  Rrt(N-l))  ♦  E(N-l) 

KHP(N)  a  EOSTROCL*  E  ( N ) ,  P  (  N ) > 

C ( N }  a  fQSTC«(L»  E  ( N ) *  RHP ! N ) *  P(N)) 
lull  fORMATUH  »«KHP! Nlaa.EiS.a*"  C  (  N  )  a"  *  E  1 5 «  b  ) 

1 1  fl  FORMAT!  lHO»"RHPa",E13P»MEa"»E13.6»"Ua",El3. 6*  "C 

la"*E13.6#"RHCa»*,E13.6*"PRHa**,E13.6) 

I F  <  AbS((RHPtN)  -  RH!N> )/RH!N> )  *GT*  TOL)  GO  TO  20 
GO  TO  22 
20  CONTINUE 

IFU  .LT.  MAX)  01)  TO  2* 

PRINT  100.  N 
GO  TU  2? 

C  GO  TO  30 

?1  CONTINUE 

IKK  *GE*  IS)  GU  TO  21* 

R  H ( N ) a (  RH(N)*RHP(N))  /  2, 

GO  10  21? 

2 1 1  FF«UP-0(N) 

FFRMaO.SEO*SlGN*C(N)*(P(N)-P(N-l ) )/RHC**2 
KH(N)aRH< N)-FF  /  FFRH 
1 F  ( UP  *lT.  0.1E-6)  GO  (0  2111 
ff»ff/up 
2  1  1  1  CONTINUE 

T  aL*TUL* 1 .E-2 

I F (  AdStFF  )  *LT.  TAD  «0  TO  22 
1010  FORMAT!  1H  » "RH ! N  )  a* ,E P « 8 # "FFa" * E 1 S ♦ 8 » "FFRH*" # E 1 5 . 8 ) 

212  continue 

C!N)aEQSTCQ(L»E!N)»RH!l',)*P!N)) 

UP«U( N  ) 

K  «  K  ♦  1 
GO  TO  10 

2?  continue 

X  (  N  )  3  xz  ♦  !  U !  N  )  ♦  SPN*C!N))*UT 
(  PRINT  900 

continue 

PF  «  P!INF) 

UF  a  U ! INF) 

110  FORMATUH  PF  a  "PlS.B*"  UF  a  "#E1*>*8) 

TO  RETURN 
LNU 
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SuBROUT I NE  REFLCKdSl.  XUU) 

C  -  SUBROUTINF  TO  CHECK  REFLECTION  OF  SECOND  SHUCK. 

COMMUN/GAIN/  0 ( 2 » 1 OOO ) » X ( 2 . 1 000 ) . UC 2 . 1 000 ) , C (  2 . 1 000 ) . KH (  •> »  1 
1000).  E(?» 10QQ)»P(2»10Q0) 

COMMON/TIMUU/  or.  UU1.  UU2»  I.XMU 
XUU  ■  X(  1.  I S 1  )  ♦  UUi*Oj 

return 

end 


subroutine  reflskid 

C  -  SUBROUTINE  FOR  SHOCK  REFLECT  AT  THE  CENTER. 

COMMON/GAIN/  0( 2. 1000) *X<2»1000) »U(2.1000) .C( 2» 1000).KH(?.  1 
1000).  E(?»1000>»P(2»1000) 

COMMON/TIMUU/  OT,  UU1.  UU2.  I.  XMU 
COHMON/NOCON/  IS1.  IS2*  I S  3 »  ISO.  INTI.  INT2.  IMAX 
COMMON  / SHK 0 A  /  U4 A . C 4 A , RH4 A . E4 A . P4 A . X4 A . GO A 
COMMON  /  SHK4H  /  U4B.C4fl.RH4B.E4B.P4B.X4B.Q4B.UUU 

common/shki/  EP 

aoo  FoRMATC"  SHOK  BEING  REFLECTED  AT  THE  CENTER”) 
rt o  1  FORMAK"  MACH  NO  FOR  fHE  REFLECTED  SHOCK  a",  ElS.fl.  /.  " 
IVELUCITy  of  the  second  shock  ■"»  Eis.b) 

rt02  FQHMATl 1H  .  12.  IX.  14*  3X»  6(E15.8»  2X).  12) 

moJ  FQRMATC1H  .  12.  IX.  I4»  3X.  6( E 1 5 . 8 »  2X).  "  REFLSHK") 

L  *  1 

XM  a  -UUl/C( 1 .  1 S2 ) 

uul  ■  "UU1 

X ( 2 .  15  2)  a  X( 1 .  I S 2  ) 

P ( 2 .  IS?)  -  P( 1 .  1S2) 

U( 2,  IS?)  a  U( 1 .  IS2) 

Rh(2»  I S2 )  *  RH ( l *  IS?) 

E  (  2 .  I S  2 )  ■  t(  1 .  IS2) 

C ( 2 .  IS?)  ■  CC 1 »  IS2) 

0(2.  IS?)  *  0(1.  IS?) 

U4A  ■  U ( 2 .  IS2) 

C 4 A  *  C ( 2 .  IS2) 

RH4A  ■  RH(  2 »  1S2) 

E4A  ■  E ( 2 .  I S 2 ) 

P4A  a  P ( 2 »  I S2  ) 

call  shqkequ.  i.  3) 

X(2*  ISi)  a  X(?»  IS2) 

P (  2 »  ISi )  r  P49 
U(2»  ISi)  a  04 B 
HH(2»  ISI)  a  RH4B 
E ( ? .  1S1 )  *  E 4fl 
C( 2.  ISi )  »  C4B 
0(2.  1S1)  a  0(2.  IS?) 

IS1H1  a  ISI  -  1 
1)0  10  K  «  1,  IS1M1 
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X(  2#  K)  *  X!  1  »  K) 

P<2»  K )  ■  P(  2 »  1 S  1  > 

HH!  2.  K)  ■  RH( 2#  iSl) 

E<  2.  K)  ■  £!  2.  ISl) 

C(2»  KJ  ■  C(  2 »  ISl) 

0(2*  K)  ■  0(2*  ISl  ) 

10  continuf 

XIS1M1  *  ISl  -  1 

UU  *  U(2»  ISl  )/XlblMl 

U{ 2  *  1  )  »  0. 

00  IS  K  *  ?,  IS1M1 
lb  U( 2*  K)  *  U(  2  »  K-l )  ♦  UU 
00  21  J  •  ISl.152 

2l  PRINT  803*I*J*X(2*J)*Ul?*J)»C(?*J)»P(2»J)»RH!2»J).E(?»J) 
HF  TURN 
EnO 


SUBROUTINE  SHFPTCL.  II*  12.  13.  14.  15.  |6.  UU) 

COMMUN/GAIN/  Ul  2.  1 000  )  *  X 1 2 » 1 000 ) , U(  2 »  1000) .  C ( 2 #  1 000 ) . RH (  ?.  1 
1000).  E(2# 10UO) »P! 2* 1000) 

COMMON  /SHK4A  /  U  4  A  .  C  4  A  ,  RH4  A  .  E  4  A  .  P4  A  »  X  4  A  »  04  A 

COMMON  /  SHK4B  /  U4H.C4fl.HH4H.E4B.P4H.X48.Q4H  .UUU 

CUMMUN/SHKI/  £  P 

CUMMON/TIMUU/  Of.  UUl.  UU2.  I.  X  MU 
commun/gnpshh/  xh 
COMMON  /  REFL  /  THEF.T 

(  dnO  FORMAT!"  XH  ■  "»  £15.8*  "  UB  ■  "»  E15.8.  «  C8  «"»  £15.0." 

C  X4H 

C  1  £15. 8) 

C  00 1  FORMAT!"  CALCULATION  Up  ITERATION  IN  SHOCK  FRONT  POINT") 

C  80 2  FORMAT!"  P«H  «"»  £15.8*  "  U4B  *"»  E15.8.  "  E4H  E15.0. 

C  1  "  RH4H  E15.8.  "  tftB  E15.8.  "  X4B  E15.8) 

C  805  FORMAT!"  SHUCK  UUES  NUT  EXlSTI  INITIATE  SHOK") 

C  8 1 0  FORMAT!"  MS  BEING  CHANGED.  MS  ■"»  12) 

020  FORMAT!  1 H  . 4X » "SHUCK  VELOCITY  FOR  THE  LEFT  SHOCK  UUl  *".£l 
15.8) 

021  FORMATUH  .  4X. "SHOCK  VELOCITY  FOR  THE  RIGHT  SHOCK  UU2  =". 
1E15.8) 

050  FORMAT! 1H  .  6II4.2X)) 

880  FORMAT  !  1 H  »  6!E15,8)J 

roo  Format  !  i  h  .  12.  lx.  u»  3x.  ME15.8.  2X).  12.  2x.  "  sfm 
IOoO  FuRMAT!"  SHUCK  FRONT  POINT  NOT  CONVERGING.") 

UINP!  V?  »  OVl.  0 V2 .  OY)  «  V 2  ♦  !  DV2*DX 1 **24DV 1 *()X?**2 ) *  OY 
1/  (  UX1*DX2*(0X1  ♦  DX2 ) )  ♦  ( "DV 1*0X2  ♦  DV2*DX 1 ) * C  OY 

2  ) * *2/  1|>X  1  *DX2  *10X1  ♦  0X2)) 

C  “  MS  ■  1,  RIGHT  RUNNING  SHOCK)  MS  ■  2.  LEFT  RUNNING  SHOCK • 

C  IF  !  UU  ,GE .  0. )  MS  «  1 
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l 


C  IF  C  UU  ,LT.  0.  >  MS  •  i 

V  IFCL  .to.  2)  00  TO  3 

C  IF(T  .  GT  *  TH£F  .AND.  *■  .EQ.  1)  GO  TO  3 

(.  1  f  ( UU  ,LE.  0.)  GO  TO  3 

C  MS  ■  2 

C  SIGN  ■  -1. 

C  Gil  TO  50 

C  3  CONTINUE 
C  Mk  ■  MS 

MS  «  2 

i  F  C  U(  1 .  13)  .LT.  Uii )  MS  ■  1 
IF(UU»  13)  .GT.  UU)  MS  •  2 
T 01.  ■  0,0005 
TuLl  *  l.E-10 
C  NIT  «  20 

N I  T  »  3  0 
tP  *  U.01 
ICHCK  *  1 
K  «  1 

c  -  uefinf  properties. 

5  GUNTINUE 


I F (  MS  . F. 0 .  1) 

SIGN  • 

1  • 

1F(MS  « F i; •  ?) 

SIGN  m 

-1. 

I F ( MS  • FO .  ?) 

GU 

TO  7 

X?  «  XU.  I?) 

U2  s  UU.  I?) 

C 2  «  CU.  I?) 

t?  «  til,  I?) 

P<?  «  Mil.  I?) 

RH2  ■  RHC 1,  12) 

02  ■  fill,  12) 

0X1  »  X( 1 .  12) 

- 

X(  1  , 

in 

DU1  ■  U( 1 .  12) 

* 

U(  1  , 

M> 

DC1  ■  C( 1 ,  12) 

m 

C(  1  , 

in 

DEI  >  Ell.  12) 

m 

E(1  . 

in 

orh i  ■  rhci,  12) 

-  RHC  1#  ID 

DPI  »  P{ 1 .  12) 

- 

PCI. 

II) 

U Q 1  *  0(1.  12) 

- 

U(  1  . 

in 

0x2  «  XU#  13) 

- 

X(  1  , 

12) 

0U2  >  U( 1.  13) 

- 

U(  1, 

1 2 ) 

DC2  «  C( 1 , 1 3  ) 

m 

C(  1, 

IO 

0E2  ■  Ell,  13) 

- 

E(l. 

1 2 ) 

0RH2  *  RH ( 1 ,  13) 

-  RHC 1,  12) 

0P2  ■  P(l,  13) 

- 

PC  1  » 

1 2 ) 

002  •  0(1#  IS) 

• 

0(1 » 

1 2 ) 

G(J  TU  « 

CONTINUE 

X2  ■  XCl,  15) 

0?  ■  UU,  15) 

C2  *  CC1,  1 5  ) 

E 
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t2  ■ 

Ed. 

15) 

P2  ■ 

P(l# 

15) 

RH2 

■  RH(  1 

•  15) 

02  ■ 

U(  1  # 

15) 

oxl 

■  X(  1# 

15) 

m 

X(  l  • 

U) 

OUl 

■  U(  1# 

15) 

m 

U(  1# 

U) 

0C1 

■  C(  1 » 

I§) 

m 

C(  1# 

U) 

0E1 

■  Ed# 

15) 

m 

Ed# 

1 4 ) 

URH1 

■  RH( 1 »  15) 

-  RH( 1# 

DPI 

■  P(  1  # 

15) 

m 

P(l» 

1  4 ) 

Ool 

■0(1# 

15) 

m 

0(1# 

1  A ) 

0X2 

■  X(  1  » 

16) 

m 

X(  1  • 

1 5 ) 

DU2 

■  U(  1 » 

16) 

m 

U(  1. 

15) 

0C2 

■  C(  1# 

16) 

m 

C(l# 

15) 

DE2 

■  E(  1 » 

16) 

m 

E(l# 

15) 

0RH2 

■  RH( 1 #16) 

-  RHd »  I 

UP2 

*  PCI* 

16) 

m 

P(  1# 

15) 

002 

■0(1# 

16) 

m 

0(1# 

15) 

C  -  COMPUTE  X4« 

8  CONTINUE 

X40  *  X( 1#  13)  ♦  UU*DT 
UB  ■  U2 
CB  ■  C 2 
Xd  *  X 2 

C  -  COMPUTE  XB  BY  ITERATION. 

9  CONTINUE 

UXB  *  xa  •  X2 

XBP  *  X 4 H  -  (UB  ♦  SlGN*CB)*DT 
Ub  ■  UlNP(U?»  DU  1  *  DU2  *  DXB> 

CB  ■  QlNP(C2#  OCl»  OC2»  DXB) 

I F  (  AdS((XBP  -  xB)/XBP>  .LT#  TOL)  GO  TO  10 
XB  ■  (XB  ♦  xBP)/2. 

C  PRINT  800.  XB.  UB.  CB#  X4B 

GO  TO  9 

10  CONTINUE 

C  I F ( MS  » F  Q  •  2)  GU  TO  20 

C  I F ( XB  .GEi  X( 1 »  13) >  G J  TO  30 

C  GO  TO  22 

(.  20  I  F<  XB  .LE.  XU#  13))  GU  TO  30 

(.  22  CONTINUE 

c  print  aoi 

U«B  ■  u? 

? 3  continue 

I F(  MS  • E0 •  1)  GO  TO  29 

CALL  GNPSH8  (L»  I  4 # 1 5 # 1 6 # X4B # MS  ) 

GO  TO  31 

?9  call  gnpshbu#  u»  12#  I3»  X4B#  ms) 

11  U48I  *  U4B 

24  IF( MS  •  FQ «  1)  GU  TO  25 

CALL  GNPSHAtL#  II#  12#  13#  X4B ) 
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GO  TO  26 

?5  call  qnpsmau*  14*  15.  16*  x«a) 

?6  call  shokeo<l»  ms*  n 

U ( (Ka(K/2)*2)  . Eu#  0)  GO  TO  265 
UUUA*UUU 
U46A«U4B 
GO  TO  266 
2b 5  UUUB-UUU 
U4BB-U4B 

2bd  IF(K  iLT » 10)  GO  TU  269 

IF( ABS(UUUA-UUUB)  .LT.u.lE-07)  GO  TO  269 
BM«( U4BA-U4HB)  /  (UUUA"UUUB) 

U4H»U4dH+  0,5*( UUUA-UUufi)*BM 
2h9  I F  (  A d S ( U4B 1 )  *LT.  TOL1)  GO  TU  27 

1 F (  A«S((U48l  •  U4B)/U4Bl>  *LT.  TOL)  GO  TO  60 

pi  continue 

IF(K  ,GE.  N I T )  GO  TO 
U4B  ■  ( U4B 1  ♦  U4B)/2. 

K  ■  K  ♦  1 

C  PRINT  602  *  P  4  6 »  U4B »  E^B*  RH4B*  C4B*  X4B 

GO  TO  23 
pa  print  iooo 
GO  TO  50 
C  30  CONTINUE 

i  print  6o5 

I F ( L  •  E  0 .  2)  GO  TO  50 
IF  l  T  ,LT.  TREE  )  MS  *  2 
IF  (  T  .GE,  TREF  )  MS  ■  1 
IF  ( MS  .EQ.  1 )  SIGN  a  1  • 

I F ( MS  .EQ.  2)  SIGN  a  -1  * 

GO  TO  50 

C  IFCL  .EO.  2)  GO  TO  50 

C  -  SMOCK  DUES  NOT  EXIST* 

(  1 F (  I C  HCk  .EO.  2)  GO  TO  50 

(  1FCMS  .EO.  1)  IS  *  13 

C  IF(MS  .EO.  2)  IS  *14 

(.  40  CONTINUE 

C  l C HCK  *  2 

C  P4A  a  P( 1*  IS) 

c  04  A  ■  U(  1  .  I  S  ) 

C  E 4 A  ■  E( I  *  IS) 

C  RH4 A  *  RM( 1 4  IS ) 

C  C  4  A  «  C(  1*  IS) 

C  U  4  A  «  0(1*  IS) 

C  I F ( MK  .EQ.  1 >  MS  ■  2 

C  IF  (  MK  . EQ .  2)  MS  a  1 

C  PRINT  610*  MS 

C  CALL  SHQKlN(L*  MS*  1*  A  S ) 

C  GO  TO  5 

so  continue 
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I F  ( L  .Fa,  2)  GU  TU  47 
uul  o  U{ 1 • I  a  > ♦  S  I  GN*C (  *.13) 

GO  TO  48 

4 7  UU2  *  U ( 1  *  I  3 ) ♦  SIGN*C(A.  13) 

4«  X40  *  X(l#  13)  ♦  (U(l*  13)  ♦  S I GN*C ( 1 »  13) ) *D T 
IF  l  L.EO.  ?)  GO  TO  55 
IF  (  MS  .Fa ,  1)  GO  TO  *5 
call  gnpshau.  ii*  l?»  13.  xat* ) 

GQ  TO  58 
55  CONTINUE 

Call  l,NPSHA(L,  14,  15.  16.  X4B) 

5«  continue 

P4d  ■  P4A 
U4B  ■  U4A 
L  4  B  «  E  4  A 
«H4B  «  RH4  A 
C  40  *  C  4  A 
040  *  Q  4  A 
ftO  LONIINUE 

LP  *  l P4H  -  P4A)/P4* 

I F  ( L  .Fa.  1  .AND.  MS  .‘■a,  2)  GO  TO  70 
X ( 2 »  13)  *  X4B 
0(2.  13)  ■  a 4 ri 
P  (  2 »  13)  *  P 4 ri 
F ( 2 .  13)  ■  F40 
KH(2»  13)  »  HH40 
C  (  2.13)  *  C  48 
U ( 2 .  13)  ■  U40 
ft  7  X ( 2  .  14)  a  X4d 

0(2.  14)  a  fl  4  A 
p ( 2  »  14)  8  p  4  A 
F ( 2.  14)  a  E4 A 

Hh(2.  la)  »  WH4  A 
C (  2  *  14)  a  C4A 
0(2.  14)  a  U4A 
GO  TU  75 
70  Continue 

X ( 2 .  13)  a  X  4  8 
0(2.  13)  ■  0  4  A 
P ( 2  »  13)  a  P4  A 
E ( 2 .  13)  «  F 4 A 
HH(2.  1 3 )  ■  HH4 A 
C ( 2 .  13)  a  C  4 A 
0(2.  13)  a  U4A 
X ( 2.  14)  a  x4d 
0(2.  14)  a  q4 A 
P ( 2 .  14)  a  p«b 
F ( 2 .  14)  a  f  40 
HH(2.  1 4 )  a  HH4d 
C ( 2 .  14)  a  C4d 
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0(2*  14)  ■  U4d 

75  continue 

IF  (  AHS(EP)  •  GT •  0.  )  GOTO  79 
I F ( Ms  » EO  *  2)  GO  TU  78 
X ( 2  #  Ii)  ■  X ( 2  *  14) 

0(2*  13)  ■  <3(2*  14) 

P ( 2  *  13)  ■  P(2»  14) 

E ( 2 *  13)  ■  E(  2*  14) 

KH(2*  13)  ■  HH( 2#  1 4 ) 

C ( 2  »  13)  ■  C ( 2 »  14) 

U ( 2 •  13)  *  U( 2*  14) 

GO  TO  79 

76  CONTINUE 

X ( 2  *  14)  «  X ( 2 •  13) 

0(2*  14)  ■  0(2*  13) 

P l 2  #  14)  *  P( 2  *  13) 

L ( 2  #  14)  *  E ( 2  *  13) 

HH(2,  14)  *  HH ( ? »  13) 

C( 2*  14)  *  C ( 2 »  13) 

U( 2*  I  4  )  *  U ( 2 »  13) 


79 

CONTINUE 

IF  (  L  . 

GO  TO  uu 

00 

Coni inue 
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CONTINUE 

IF  (  L  .EQ 
IF  (  L  » E  U  • 

return 

end 


1  )  GO  TO  do 


1  )  uu  »  Pul 

2  )  UU  «  uU? 


subroutine  shukin  (L*M&«Mg* is)  ^ 

common/ gain/  q(2*iooo)*x(2»iooo)»u(2» iooo )»ci2*iooo)#rh(?*i 
1000)*  EC2»1000)#P(2»1000) 

COHMUN/NOCQN/  ISl*  IS?*  I S  3 •  IS4*  INTI#  INT2*  IMAX 
COMMON  /  SHK4H  /  U48.C4B»HH4B*E4B.P4B*X4B*Q48  »UUU 

Common  /SHK4A  /  U4A»C4A#RH4AtE4A*P4A*X4A»Q4A 
C0MM0N/TI MUU/  0T»  UU 1  *  UU2»  I*  XMU 
cqmmun/shki/  ep 

c  is  position  to  initiate  smock 

c  MS  a  1  -  RIGHT  SHOCK*  MS  *  2  "  LEFT  SHOCK 

C  mQ  *  "  GIVEN  P  2  *  MQ  a  2  -  GIVEN  U2 

1000  KQRMATUH  » I  2*  IX*  14. 3X*  ft(E15«8»2X) » I2*2X»"INS«"  ) 

100 1  Format ( *  initiated  shuck  speed  uu2  ■"»  E15.8) 

1002  format ( w  initiated  shuck  speed  uui  ■  "»  eis.8> 

ToESK  *  1*E»06 

Tolsk-o.oooi 

IF  (  MS  .EQ.  1)  SIGN  *  1* 

IF  (  MS  .EQ.  2  )  SIGN  *  “1 « 
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151  «  IS 

1 5 2  ■  Is  ♦  1 

U4A  *  U( 1  *  IS) 

C4A  *  C(  1*  IS) 

NH4A  *  RHC  1 »  IS) 

UA  ■  E<  1.  IS) 

P4A  ■  P<  t,  IS) 

04 A  ■  0(1*  IS) 

IF  C  AUS(EP)  .LT.  TOLS*  )  GO  TO  7 
P4B  *  P4*  *  EP  ♦  PA A 
CALL  SHOKEQ  CL. MS*  1) 

(iO  TO  a 
7  CONTINUE 

Uul  ■  U(l.IS)  ♦  SIGN  *  C (  1  *  I S ) 
CALL  SHlIKEO  (  L.MS.  3) 
b  CONTINUE 

IFCMS  .EG.  1)  GU  TO  10 


U(l.lsl) 

U4  A 

PC  1 # IS1) 

P«A 

RH  Cl.ISl 

■  RH4  A 

E(I.ISI) 

E4A 

C( 1 , IS1  ) 

C4A 

0 ( 1.IS1) 

04  A 

X( 1 .  IS1  ) 

XC1.IS) 

U( 1 .  IS2) 

U4ti 

P( 1.IS2) 

P4U 

HH( 1 . IS?) 

■  HH4b 

t( 1.IS2) 

E4b 

C( 1 .  IS2) 

C4d 

U( 1 . IS2) 

OC 1  *  IS) 

X( 1  * IS2) 
GO  TO  20 
CONTINUE 

XC 1 .  IS) 

U( 1 » IS1 ) 

U4B 

P(l.ISl) 

P4d 

HHC  1.  IS1 ) 

■  HH4d 

E  (  1  *  I  S 1  ) 

E  4B 

C( 1.IS1) 

C4« 

Q( 1 # ISl  ) 

0C1.1S) 

X( 1  .  IS1  ) 

XC1.1S) 

U( 1  *  IS2) 

U4A 

P(  1.IS2) 

P4  A 

HHC  1. IS?) 

*  HH4A 

f(  1.IS2) 

f  4  A 

C( 1* IS2) 

C4  A 

OC 1# IS2) 

OC  1  *  IS) 

X( 1 . IS2) 
CONTINUE 

XC 1  *  IS  ) 

1ECL  ,E0, 
GO  TO  2? 

1)  GO  To 
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21  continue 
? i  continue 
return 
end 


SuQROuTINE  SHOkEQ  (I.mS.MQ) 

COMMON  /SHK4A  /  U4A.C4A.RH4A.E4A.P4A.X4A.04A 
Common  /  SHK4B  /  U4B.C4fl.RH4B.E4B.P48.X4B.Q4B.UUU 
COMMON/TIMUU/  OT.  UUl.  UU2  *  I#  XMU 
COMMON/SHKI/  ep 
C  9()U  F  URMAT  (  lrt  »  6(  £15. ft)  > 

C  MOl  F  OHM  A  T  (  1 H  ROUTINE  SMOKEO  L«">13»"  MS*". 12."  MO*". 14) 

f.  MS  *  1  RIGHT  SHO^K 

C  MS*  2  -  LEFT  SHOCK 

mu  *  1,  GIVEN  P2J  MU  *  ?.  GIVEN  1)2)  MO  *  J.  GIVEN  DO, 

MU  *  1  -  GIVEN  PP»  Mg  *2  -  GIVEN  U2 

PRINT  410 

C  410  FORMAT  C1H  ,  "  U2P  UUP  E2P  RH2P  C2P  HP") 

406  FUHMAl  (1H  .  "  ***  WARNING  ****  RH2  IS  LESS  THAN  RMi"  ) 

looo  format  c i h  ."shuckeo  nut  converging") 

C  PRINT  901.L.MS.MQ 

1 j  1  3  f 0RMAT(6X."RMl»  ".E15. °*6X»"UU*  "  » E  1 5 » 8 » 6X  # "  U 1  *  ".E15.8) 

xs  *  1  . 

IF  (  MS  • EO .  2)  XS  *  -I. 

TQL  *  0,0005 
T0L1*  1  .E-20 
N  I  T  *  1  0  0 
K*  1 

U 1  *  U  4  A 
C I*C4 A 

uu  ■  Cl  *  xs  ♦  u  1 

HHl*RM«A 
t  l*t4A 
P  1  *  P  4  A 

GQ  TU  (5.  25.  40).  MQ 
5  continue 
C  PRINTV03 

C  403  FORMATUH  .6X."U2P",12A»"UUP"»12X»"E2P"#12X»"KH2P".11X» 

C  "C2P".12X. 

C  *"P2P") 

P2*P4d 

p?p  a  pan 

I F l ( P2  -  Pi)  .LE,  0.)  P?  ■  P 1  ♦  0.001 
U2  *  Ul  ♦  (P2-P1 )/(RHl*(UU-Ul )  ) 

RH2*C  P2/P1 )**1 ,4*HH1 
10  CONTINUE 

IF  (  A8S(U2-U1)  .LT.  T^lI)  GO  TO  15 

U2P  «  (P2-PI)*(RH2“RH1}/(RH1*HH2*(U2-U1))+U1 
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.  ***•  s  « ^&m**#*'*  ***********  «— 


12  UUP«(P2-Pl)/(KHl«(U2P-UimUl 

15  IF( (RH2-RH1 ) »LE*0.0)CALl  EXIT 

E2P  «  El  ♦  <PDP2)*(RH*-RH1  ) / < 2  * Q*RH 1 *RH2 ) 
HH2Pafc.«sTR0(L#E2P#P2) 

C2P«EUSTCU( L#E2P#RH2P»P2) 

IF  (  ABS  ((U2P  -  U2  )  I  U2P)  -TOL  )  16#  16.  20 

16  IF  (ABS<(  RH2P  “  RH2  )'RH2P)  -  TOL  )  30#30,20 
?0  IF  (  K  ,GE.  NIT  )  GO  TO  35 

I F ( (  K-(K/2)*2)  .EO.  0)  GO  TO  202 
RHX A*RH? 

RHYA-RH2P 
GU  TU  205 
202  RHXB*Rh2 
RHYB*Rp?P 

2()5  IF  (k  •  LT  « 1 2 )  GO  TO  2^6 

IF  ( AHS(RHYB-RHYA)»LT*  0.1E-06)  GO  TU  ?0B 
BOb*  (RHXfl  -  RHXA)/(R'1f8-RHYA> 

RH?*  (RHXA-BUfl*RHYA)/I U  -BOB) 

GO  TU  208 

2(i6  HH2*0.5*(Rh2  ♦  RH2P ) 

2<)«  K»K*1 

U?*0 , 5* ( U2+U2P ) 

C  PRINT  R0O#U?P»UUP#E2P#rtH2P»C2P»P2P 

GO  TO  10 
35  PRINT  1000 
GO  TO  30 
?5  continue 

U2  b  U4B 
U?P  *  U4H 
R  H  2  *  RHl  *  1.1 
26  CONTINUE 

IF  (  ABS  (Rh2  -RHl)  .1*.  T0L1)  GO  TO  27 
UUP  •  ( RH2*U2  -  RH1*U1 )/(RH2-RHl) 

?7  P?P  ■  Pi  ♦  RH 1  * ( UUP  -  U|)*(U2aUl) 

EpP  *  El  ♦  (Pit  P2P)  *  (RH2-RH1  )/l2.*HHl*RH2) 
HH^r  b  FOSTRQ  (L»E?P*P2p) 

CpP  *  E.QSTCO  (L#E2P#RH<!p,P2P) 

IF  (  ABS((RH2P  -RH2)/RH2P)  .LT.  TOL)  GO  TO  30 
K  *  Kt- 1 

RH2  *  IRH2P  ♦  RH2)/2, 

C  PRINT  900.U?P#UUP#E2P.rtH2P#C?P»P?P 

IF  (  K  .  Gt .  NIT)  GO  TO  29 
GU  TU  26 

?v  print  1000 

G()  TO  30 
40  CONTINUE 


UUP 

«  uui 

RH? 

*  PHI  *2.0 

I  F  ( 

Ab5(EP)  .LT. 

0,001)  RH2 

■  RHl 

4  I  U?P 

a  (RH2*UUP  - 

RH 1  * ( uyP  - 

U1  )  )/RH2 
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H2P  ■  Pi  ♦  RW1*(UUP  •  ^1)*(U?P  *  Ul) 

E2P  *  El  ♦  (PU  P2P)  *  {  RH2-RH1  )  /  (  2 .  *HHl  *KM2 ) 
hh2p  ■  eqstru  <l.E2p.p*p) 

C2P  »  EoSTCo  (L#£2P.RH2p»P2P) 

IF  (  ABS((WH2P  •RH?’)/RH2P)  »LT»  TOL)  GO  TO  30 


K  ■  K+l 

PRINT  900»U?P#UUP#E2P.«H2P.C?P#P2P 
IF  l  K  .G£.  NIT)  GO  TO  39 
KH2  «  (RH2P  ♦  RM2)/2. 

GO  TO  41 
print  iooo 
)0  CONTINUF 
U4B*U2P 
C4B«C2P 
RH4H«nH2P 
E  4  B  ■  E  2  P 
H 4B  *  P2P 

1^9  F[)RMAT("  UUP»".E15.8) 

IF  C  L  .Ed.  1)  UU1  «  U^P 
IF  l  L  .Ed.  2)  UU2  •  U^P 


(' 


U4B.UUP»E4H*RH4B»C4R»P«B 


uuu  *  uup 

PRINT  900. 

RETURN 

END 

SUBROUTINE  SRlTCH(II.I^On.IHAT.ISKOS)  .. 

CnMMUN/GAlN/«(?.  1000),M  2*  1000).U(2.1000).C(?.1000),RH(2.10 


100).  E(2.1000).P(2.1U00) 

COMMUN/NOCDN/ISl. 1S2.IS3. IS4.INT1 . INT2.MXNPT 


I F  (  I AUO.LT .0)  GU  TO  1 
I  S  I  GN*  1 


N«0 


GU  TO  2 

1  I S 1 GN*“  1 

n*mxnpt-i  I 

2  CONTINUE 
ilast*mxnpt*iadu-n-i 
00  3  J*II.MXNPT 

iach»mxnptmado*(  II-J)*ISI gn-n 

lQU»MXNPT>(Ii“J)*ISIGN“N 
X(1»IACH)*X( l.IUU) 

U( 1. IACH)«U( l.IOU) 

C ( 1  *  I ACH ) *C  (l.IOU) 

P( l, IACH)*P( 1 » IOU) 

RH(l»l^CH)*RH(l#IUU) 
E(1#IACH)»E( l.IUU) 
d( 1 . I ACH ) *Q( 1 # IOU) 
j  CONTINUE 

MxNPT«MxNPT*I ado 

return 

ENU 
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